A Flavor Kit for BSM models

We present a new kit for the study of flavor observables in models beyond the standard model. The setup is based on the public codes SARAH and SPheno and allows for an easy implementation of new observables. The Wilson coefficients of the corresponding operators in the effective Lagrangian are computed by SPheno modules written by SARAH. New operators can also be added by the user in a modular way. For this purpose a handy Mathematica package called PreSARAH has been developed. This uses FeynArts and FormCalc to derive the generic form factors at tree- and 1-loop levels and to generate the necessary input files for SARAH. This framework has been used to implement BR($\ell_\alpha\to \ell_\beta \gamma$), BR($\ell_\alpha \to 3\, \ell_\beta$), CR($\mu-e,A$), BR($\tau \to P \, \ell$), BR($h\to \ell_\alpha \ell_\beta$), BR($Z\to \ell_\alpha \ell_\beta$), BR($B_{s,d}^0 \to \ell \bar{\ell}$), BR($\bar B \to X_s\gamma$), BR($\bar B \to X_s \ell \bar{\ell}$), BR($\bar B \to X_{d,s} \nu \bar \nu$), BR($K^+ \to \pi^+ \nu \bar \nu$), BR($K_L \to \pi^0 \nu \bar \nu$), $\Delta M_{B_s,B_d}$, $\Delta M_K$, $\varepsilon_K$, BR($B \to K \mu \bar{\mu}$), BR($B\to \ell \nu$), BR($D_s \to \ell \nu$) and BR($K \to \ell \nu$) in SARAH. Predictions for these observables can now be obtained in a wide range of SUSY and non-SUSY models. Finally, the user can use the same approach to easily compute additional observables.


Introduction
With the exploration of the terascale, particle physics has entered a new era. On the one hand, the discovery of a Higgs boson at the LHC [1,2] seemingly completed the Standard Model (SM) of particle physics, even though there is still quite some room for deviations from the SM predictions. The observed mass of about 125 GeV in combination with a top quark mass of 173.34 GeV [3] implies within the SM that we potentially live in a metastable vacuum [4]. This, together with other observations, like the dark matter relic density or the unification of gauge forces, indicates that there is physics beyond the SM (BSM). Although no sign of new physics has been found so far at the LHC, colliders are not the only places where one can search for new physics. Low energy experiments focused on flavor observables can also play a major role in this regard, since new particles leave their traces via quantum effects in flavor violating processes such as b → sγ, Bs → µ + µ − or µ → eγ.
In the last few years there has been a tremendous progress in this field, both on the experimental as well as on the theoretical side. In particular, observables from the Kaon-and B-meson sectors, rare lepton decays and electric dipole moments have put stringent bounds on new flavor mixing parameters and/or additional phases in models beyond the SM.
There are several public tools on the market which predict the rates of several flavor observables: superiso [5][6][7], SUSY_Flavor [8,9], NMSSM-Tools [10], MicrOmegas [11][12][13][14][15], SuperBSG [16], SupeLFV [17], SuseFlav [18], IsaJet with IsaTools [19][20][21][22][23][24] or SPheno [25,26]. However, all of these codes have in common that they are only valid in the Two-Higgs-doublet model or in the MSSM or simple extensions of it (NMSSM, bilinear Rparity violation). In addition, none of these tools can be easily extended by the user to calculate additional observables. This has made flavor studies beyond the SM a cumbersome task. The situation has changed with the development of SARAH [27][28][29][30][31]. This Mathematica package can be used to generate modules for SPheno, which then can calculate flavor observables at the 1-loop level in a wide range of supersymmetric and nonsupersymmetric models [32][33][34]. However, so far all the information about the underlying Wilson coefficients 1 for the operators triggering the flavor violation as well as the calculation of the flavor observables had been hardcoded in SARAH. Therefore, it was also very difficult for the user to extend the list of calculated observables. The implementation of new operators was even more difficult.
We present a new kit for the study of flavor observables beyond the standard model. In contrast to previous flavor codes, FlavorKit is not restricted to a single model, but can be used to obtain predictions for flavor observables in a wide range of models (SUSY and non-SUSY). FlavorKit can be used in two different ways. The basic usage of FlavorKit allows for the computation of a large number of lepton and quark flavor observables, using generic analytical expressions for the Wilson coefficients of the relevant operators. The setup is based on the public codes SARAH and SPheno, and thus allows for the analytical and numerical computation of the observables in the model defined by the user. If necessary, the user can also go beyond the basic usage and define his own operators and/or observables. For this purpose, a Mathematica package called PreSARAH has been developed. This tool uses FeynArts/FormCalc [35][36][37][38][39][40] to compute generic expressions for the required Wilson coefficients at the tree-and 1-loop levels. Similarly, the user can easily implement new observables. With all these tools properly combined, the user can obtain analytical and numerical results for the observables of his interest in the model of his choice. To calculate new flavor observables with SPheno for a given model the user only needs the definition of the operators and the corresponding expressions for the observables as well as the model file for SARAH. All necessary calculations are done automatically. We have used this setup to implement BR( α → β γ), BR( α → 3 β ), CR(µ − e, A), BR(τ → P ), BR(h → α β ), BR(Z → α β ), BR(B 0 s,d → ¯ ), BR(B → Xsγ), BR(B → Xs ¯ ), BR(B → X d,s νν), BR(K + → π + νν), BR(K L → π 0 νν), ∆M Bs,B d , ∆M K , ε K , BR(B → Kµμ), BR(B → ν), BR(Ds → ν) and BR(K → ν) in SARAH.
This manual is structured as follows: in the next section we give a brief introduction into the calculation of flavor observables focusing on the main steps that one has to follow. Then we present FlavorKit, our setup to combine FeynArts/FormCalc, SPheno and SARAH in Section 3. In Section 4 we explain how new observables can be added and in Section 5 how the list of operators can be extended by the user. A comparison between FlavorKit and the other public codes is presented in Section 6 taking the MSSM as an example before we conclude in Section 7. The appendix contains information about the existing operators and how they have been combined to compute the different flavor observables.

General strategy: calculation of flavor observables in a nutshell
Once we have chosen a BSM model 2 , our general strategy for the computation of a flavor observable follows these steps: -Step 1: We first consider an effective Lagrangian that includes the operators relevant for the flavor observable of our interest, This Lagrangian consists of a list of (usually) higher-dimensional operators O i . The Wilson coefficients C i can be induced either at tree or at higher loop levels and include both the SM and the BSM contributions ). They encode the physics of our model. The user has to make a choice in step 1. The list of operators in the effective Lagrangian can be restricted to the most relevant ones or include additional operators beyond the leading contribution, depending on the required level of precision. Usually, the complete set of renormalizable operators contributing to the observable of interest is considered, although in some well motivated cases one may decide to concentrate on a smaller subset of operators. This freedom is not present in step 2. Once the list of operators has been arranged, the computation of the corresponding C i coefficients follows from the consideration of all topologies (penguin diagrams, box diagrams, . . . ) leading to the O i operators. This is the most complicated and model dependent step, since it demands a full knowledge of all masses and vertices in the model under study. Furthermore, it may be necessary to compute the coefficients at an energy scale and then obtain, by means of their renormalization group running, their values at a different scale. Finally, step 3 is usually quite straightforward since, like step 1, is model independent. In fact, the literature contains general expressions for most flavor observables, thus facilitating the final step. However, one should be aware that the formulas given in the literature assume that certain operators contribute only sub-dominantly and, thus, omit the corresponding contributions. This is in general justified for the SM but not in a general BSM model. In particular, this is the case for processes involving external neutrinos, which are often assumed to be purely left-handed, making the operators associated to their right-handed components to be neglected.
We will exemplify our strategy using a simple example: BR(µ → eγ) in the Standard Model extended by right-handed neutrinos and Dirac neutrino masses. The starting point is, as explained above, to choose the relevant operators. In this case, it is well known that only dipole interactions can contribute to to the radiative decay α → β γ at leading order 4 . Therefore, the relevant operators are contained in the − − γ dipole interaction Lagrangian. This is in general given by L dipole γ = ie m α¯ β σ µν qν K L 2 P L + K R 2 P R αAµ + h.c. (2) Here e is the electric charge, q the photon momentum, P L,R = 1 2 (1 ∓ γ 5 ) are the usual chirality projectors and α,β denote the lepton flavors. This concludes step 1.
The information about the underlying model is encoded in the coefficients K L,R 2 . In the next step, these coefficients have to be calculated by summing up all Feynman diagrams contributing at a given loop level. Expressions for these coefficients for many different models are available in the literature. In the SM only neutrino loops contribute and one finds [41] The current version of FlavorKit can only handle renormalizable operators at this stage of the computation. 3 In principle, one can go beyond the 1-loop level, although in our case we will restrict our computation to the addition of a few NLO corrections. 4 At next to leading order, one would also have to consider operators likeμγν eqγ ν q, to be combined with a q − q − γ dipole interaction.
Here, λ ij denote the entries of the Pontecorvo-Maki-Nakagawa-Sakata matrix and F 1 and F 2 are loop functions.
One finds approximately F 1 − 1 4 mν m W 2 and F 2 0. Finally, we just need to proceed to the last step, the computation of the observable. After computing the Wilson coefficients K L,R 2 it is easy to relate them to BR(µ → eγ) by using [42] Γ α → β γ = αm 5 This expression holds for all models. With this final step, the computation concludes.
As we have seen, the main task to get a prediction for BR(µ → eγ) in a new model is to calculate K L,R 2 . However, this demands the knowledge of all masses and vertices involved. Moreover, in most cases a numerical evaluation of the resulting loop integrals is also welcome. Therefore, even for a simple process like µ → eγ, a computation from scratch in a new model can be a hard work. In order to solve this practical problem, we are going to present here a fully automatized way to calculate a wide range of flavor observables for several classes of models.

FlavorKit: usage and goals
As we have seen, the calculation of flavor observables in a specific model is a very demanding task. A detailed knowledge about the model is required, including 1. expressions for all involved masses and vertices 2. optionally, renormalization group equations to get the running parameters at the considered scale 3. expressions to calculate the operators 4. formulae to obtain the observables from the operators Nearly all codes devoted to flavor physics have those pieces hardcoded, and they are only valid for a few specific models 5 . The only exception is SPheno, thanks to its extendability with new modules for additional models. These modules are generated by the Mathematica package SARAH and provide all necessary information about the calculation of the (loop corrected) mass spectrum, the vertices and the 2-loop RGEs. These expressions, derived from fundamental principles for any (renormalizable) model, contain all the information required for the computation of flavor observables. In fact, SARAH also provides Fortran code for a set of flavor observables. For this output, generic expressions of the necessary Wilson coefficients have been included. These are matched to the model chosen by the user and related to the observables by the standard formulae available in the literature. However, it was hardly possible for the user to extend the list of observables or operators included in SARAH without a profound knowledge of either the corresponding Mathematica or Fortran code.
We present a new setup to fill this gap in SARAH: FlavorKit. As discussed in Sec. 2, the critical step in the computation of a flavor observable is the derivation of analytical expressions for the Wilson coefficients of the relevant operators. This step, being model dependent, requires information about the model spectrum and interactions. However, generic expressions can be derived, later to be matched to the specific spectrum and interaction Lagrangian of a given model. For this purpose, we have created a new Mathematica package called PreSARAH. This package uses the power of FeynArts and FormCalc to calculate generic 1-loop amplitudes, to extract the coefficients of the demanded operators, to translate them into the syntax needed for SARAH and to write the necessary wrapper code. PreSARAH works for any 4-fermion or 2-fermion-1-boson operators and will be extended in the future to include other kinds of operators. The current version already contains a long list of fully implemented operators (see Appendix B). The results for the Wilson coefficients obtained with PreSARAH are then interpreted by SARAH, which adapts the generic expressions to the specific details of the model chosen by the user and uses snippets of Fortran code to calculate flavor observables from the resulting Wilson coefficients. As for the operators, there is a long list of observables already implemented (see Appendices C.1 and C.2). Finally, SARAH can be used to obtain analytical output in L A T E X format or to create Fortran modules for SPheno, thus making possible numerical studies.
FlavorKit can be used in two ways: Lepton flavor Quark flavor Table 1 List of flavor violating processes and observables which have been already implemented in FlavorKit. To the left, observables related to lepton flavor, whereas to the right observables associated to quark flavor. See appendicesC.1 andC.2 for the definition of the observables and the relevant references for their calculation.
-Basic usage: This is the approach to be followed by the user who does not need any operator nor observable beyond what is already implemented in FlavorKit. In this case, FlavorKit reduces to the standard SARAH package. The user can use SARAH to obtain analytical results for the flavor observables and, if he wants to make numerical studies, to produce Fortran modules for SPheno. For the list of implemented operators we refer to AppendixB, whereas the list of implemented observables is given in Table 1. -Advanced usage: This is the approach to be followed by the user who needs an operator or an observable not included in FlavorKit. In case the user is interested in an operator that is not implemented in FlavorKit, he can define his own operators and get analytical results for their coefficients using PreSARAH.
Then the output can be passed to SARAH in order to continue with the basic usage. In case the user is interested in an observable that is not implemented in FlavorKit, this can be easily implemented by the addition of a Fortran file, with a few lines of code relating the observable to the operators in FlavorKit (implemented by default or added by the user). The Fortran files just have to be put together with a short steering file into a specific directory located in the main SARAH directory. Then one can continue with the basic usage. The combination of PreSARAH together with SARAH and SPheno allows for a modular and precise calculation of flavor observables in a wide range of particles physics models. We have summarized the setup in Fig. 1: the user provides as input SARAH model files for his favorite models or takes one of the models which are already implemented in SARAH (see Appendix D for a list of models available in SARAH). New observables are implemented by providing the necessary Fortran code to SARAH while new operators can be either implemented by hand or by using PreSARAH which then calls FeynArts and FormCalc for the calculation of the necessary diagrams. However, most users will not require to implement new operators or observables. In this case, the user can simply use SARAH in the standard way and (1) derive analytical results for the Wilson coefficients and observables, and (2) generate Fortran modules for SPheno in order to run numerical analysis.

Download and installation
FlavorKit involves several public codes. We proceed to describe how to download and install them. In case the user is interested in the numerical evaluation of the flavor observables, a SPheno module must be created as explained above. Once this is done, the resulting Fortran code can be used for the numerical analysis of the model. This can be achieved in the following way: In those two blocks, the results for quark and lepton flavor violating observables are given.
Finally, an even easier way to implement new models in SARAH is the butler script provided with the SUSY Toolbox [45] sarah.hepforge.org/Toolbox/

Limitations
FlavorKit is a tool intended to be as general as possible. For this reason, there are some limitations compared to codes which perform specific calculations in a specific model. Here we list the main limitations of FlavorKit: -Chiral resummation is not included because of its large model dependence, see e.g. [46] and references therein.
-Even though we have included some of the higher order corrections for the SM part of some observables in a parametric way, 2-or higher loop corrections, calculated in the context of the SM or the MSSM for specific observables, are not considered, see for instance [47][48][49][50][51][52][53][54]. The steering file contains the following information: -NameProcess: a string as name for the set of observables.
-NameObservables: names for the individual observables and numbers which are used to identify them later in the SPheno output. The value is a three dimensional list. The first part of each entry has to be a symbol, the second one an integer and the third one a comment to be printed in the SPheno output file ({{name1,number1,comment1},...}). -NeededOperators: The operators which are needed to calculate the observables. A list with all operators already implemented in FlavorKit is given in Appendix B. In case the user needs additional operators, this is explained in Sec. 5. -Body: The name (as string) of the file which contains the Fortran code to calculate the observables from the operators.
For instance, the corresponding file to calculate α → β γ reads The observables will be saved in the variables muEgamma, tauEgamma, tauMuGamma and will show up in the spectrum file written by SPheno in the block FlavorKitLFV as numbers 701 to 703.
Besides the operators, the SM parameters given in Table 2 and the hadronic parameters given in Tables 3  and 4 can be used in the calculations. For instance, we used Alpha for α(0) and mf_l which contains the poles masses of the leptons as well as GammaMu and GammaTau for the total widths of µ and τ leptons.

Real Variables
AlphaS_MZ Table 2 List of SM parameters available in FlavorKit. All hadronic observables are calculated at Q = 160 GeV.
By extending or changing the file hadronic_parameters.m in the FlavorKit directory, it is possible to add new variables for the mass or life time of mesons. These variables are available globally in the resulting SPheno code. The numerical values for the hadronic parameters can be changed in the Les Houches input file by using the blocks FCONST and FMASS defined in the Flavor Les Houches Accord (FLHA) [55].
It may happen that the calculation of a specific observable has to be adjusted for each model. This is for instance the case when (1) the calculation requires the knowledge of the number of generations of fields, (2) the mass or decay width of a particle, calculated by SPheno, is needed as input, or (3) a rotation matrix of a specific field enters the analytical expressions for the observable. For these situations, a special syntax has been created. It is possible to start a line with @ in the Fortran file. This line will then be parsed by SARAH, and Mathematica commands, as well as SARAH specific commands, can be used. We made use of this functionality in the implementation of h → α β .  returns the name used for the mass of a particle x in the SPheno output SPhenoMassSq [x] returns the name used for the mass squared of a particle x in the SPheno output SPhenoWidth [x] returns the name used for the width of a particle x in the SPheno output HiggsMixingMatrix name of the mixing matrix for the CP even Higgs states in a given model PseudoScalarMixingMatrix name of the mixing matrix for the CP odd Higgs states in a given model Table 5 SARAH commands which can be used in the input file for the calculation of an observable. We give in Table 5 the most important SARAH commands which might be useful in this context.
Many more examples are given in AppendixC.1, where we have added all input files for the calculations of flavor observables delivered with SARAH.

Advanced usage II: Implementation of new operators
The user can also implement new operators and obtain analytical expressions for their Wilson coefficients. In this case, he will need to use PreSARAH which, with the help of FeynArts and FormCalc, provides generic expressions for the coefficients, later to be adapted to specific models with SARAH.

Introduction
New operators can be implemented by extending the content of the folder

[$SARAH]/FlavorKit/[$Type]/Operators/
In the current version of FlavorKit, 3-and 4-point operators are supported. Each operator is defined by a .m-file. These files contain information about the external particles, the kind of considered diagrams (tree-level, self-energies, penguins, boxes) as well as generic expressions for the coefficients. These expressions, derived from the generic Feynman diagrams contributing to the coefficients, are written in the form of a Mathematica code, which can be used to generate Fortran code.
For the automatization of the underlying calculations we have created an additional Mathematica package called PreSARAH, which can be used to create the files for all 4-fermion as well as 2-fermion-1-boson operators. This package creates not only the infrastructure to include the operators in the SPheno output of SARAH but makes also use of FeynArts and FormCalc to calculate the amplitudes and to extract the coefficient of the demanded operators. It takes into account all topologies depicted in Figs. 2 to 6.

Input for PreSARAH
In order to derive the results for the Wilson coefficients, PreSARAH needs an input file with the following information: is done by using the Description statements Higgs or Pseudo-Scalar Higgs as well as Scalar-Mixing-Matrix or Pseudo-Scalar-Mixing-Matrix. If the particle or parameter needed to calculate an observable is not present or has not been defined, the observable is skipped in the SPheno output.  -ConsideredProcess: A string which defines the generic type for the process -"4Fermion" -"2Fermion1Scalar" -"2Fermion1Vector" -NameProcess: A string to uniquely define the process -ExternalFields: The external fields. Possible names are ChargedLepton, Neutrino, DownQuark, UpQuark, ScalarHiggs, PseudoScalar, Zboson, Wboson 7 -FermionOrderExternal: the fermion order to apply the Fierz transformation (see the FormCalc manual for more details) -NeglectMasses: which external masses can be neglected (a list of integers counting the external fields) -ColorFlow: defines the color flow in the case of four quark operators. To contract the colors of external fields, ColorDelta is used, i.e ColorFlow = ColorDelta[1,2]*ColorDelta[3,4] assigns (q α Γ qα)(q β Γ q β ).
-AllOperators: a list with the definition of the operators. This is a two dimensional list, where the first entry defines the name of the operator and the second one the Lorentz structure. The operators are expressed in the chiral basis and the syntax for Dirac chains in FormCalc is used: for the helicity of an external gauge boson.
k[N] for the momentum of the external particle N (N is an integer). 7 The particles.m file is used to define for each model which particle corresponds to SM states using the Description statement together with Leptons, Neutrinos, Down-Quarks, Up-Quarks, Higgs, Pseudo-Scalar Higgs, Z-Boson, W-Boson.
If there is a mixture between the SM particles and other states (like in R-parity violating SUSY or in models with additional vector quarks/leptons) the combined state has to be labeled according to the description for the SM state.
Notice that in the SM Pseudo-Scalar Higgs is just the neutral Goldstone boson. If an external state is not present in a given model or has not been defined as such in the particles.m file the corresponding Wilson coefficients are not calculated by SPheno.  -A Dirac chain starting with a negative first entry is taken to be anti-symmetrized. See the FormCalc manual for more details. To make the definitions more readable, not the full DiracChain object of FeynArts/FormCalc has to be defined: PreSARAH puts everything with the head Op into a Dirac chain using the defined fermion order. For 4-fermion operators the combination of both operators is written as dot product. For instance Op [6].Op [6] is internally translated into   { OllddSLL , Op [ 7 ] . Op [ 7 ] } , 13 {OllddSRR , Op [ 6 ] . Op [ 6 ] } , 14 {OllddSRL , Op [ 6 ] . Op [ 7 ] } , 15 {OllddSLR , Op [ 7 ] . Op [ 6 ]   Here, we neglect all external masses in the operators (NeglectMasses={1,2,3,4}), and the different coefficients of the scalar operators (¯ P X )(dP Y d) are called OllddSXY, the ones for the vector operators (¯ P X γµ )(dP Y γ µ d) are called OllddVYX and the ones for the tensor operators (¯ P X σµν )(dσ µν P Y d) OllddTYX, with X,Y=L,R.
Notice that FormCalc returns the results in form of P X γµ while in the literature the order γµP X is often used. Finally, SPheno will not calculate all possible combinations of external states, but only some specific cases: µedd, τ edd, τ µdd, µess, τ ess, τ µss 8 .
The input file to calculate the coefficients of the − − Z operators (¯ γµP L,R )Z µ and (¯ pµP L,R γµ )Z µ is 1 NameProcess="Z 2 l " ;   In this case the operators for all Higgs states present in the considered model will be computed.

Operators with massless gauge bosons
We have to add a few more remarks concerning 2-fermion-1-boson operators with massless gauge bosons since those are treated in a special way. It is common for these operators to include terms in the amplitude which are proportional to the external masses. Therefore, if one proceeds in the usual way and neglects the external momenta, some inconsistencies would be obtained. For this reason, a special treatment is in order. In PreSARAH, when one uses the dependence on the two fermion masses is neglected in the resulting Passarino-Veltman integrals but terms proportional to m f1 and m f2 are kept. This solves the aforementioned potential inconsistencies. Furthermore, for the dipole operators, defined by we are using the results obtained by FeynArts and FormCalc and have implemented all special cases for the involved loop integrals (C 0 , C 00 , C 1 , C 2 , C 11 , C 12 , C 22 ) with identical or vanishing internal masses in SPheno. This guarantees the numerical stability of the results 9 .
The monopole operators of the form q 2 (f γµf )V µ are only non-zero for off-shell external gauge bosons, while PreSARAH always treats all fields as on-shell. Because of this, and to stabilize the numerical evaluation later on, these operators are treated differently to all other operators: the coefficients are not calculated by FeynArts and FormCalc but instead we have included the generic expressions in PreSARAH using a special set of loop functions in SPheno. In these loop functions the resulting Passarino-Veltman integrals are already combined, leading to well-known expressions in the literature, see [42,56]. They have been cross-checked with the package Peng4BSM@LO [43]. To get the coefficients for the monopole operators, these have to be given always in the form in the input of PreSARAH.

Combination and normalization of operators
The user can define new operators as combination of existing operators. For this purpose wrapper files containing the definition of the operators can be included in the FlavorKit directories. These files have to begin with ProcessWrapper = True;. This function is also used by PreSARAH in the case of 4-fermion operators: for these operators the contributions stemming from tree-level, box-and penguin-diagrams are saved separately and summed up at the end. Thus, the wrapper code for the 4-lepton operators written by PreSARAH reads  It is also possible to use these wrapper files to change the normalization of the operators. We have made use of this functionality for the operators with external photons and gluons to match the standard definition used in literature: it is common to write these operators as e m f (f σµν f )F µν , i.e. with the electric coupling (or strong coupling for gluon operators) and fermion mass factored out. This is done with the files Photon_wrapper.m and Gluon_wrapper.m, included in the FlavorKit directory of SARAH:  First, the coefficients OA2qSL and OA2qSR derived with PreSARAH are matched to the new coefficients CC7 and CC7p. The same matching is automatically applied also to the SM coefficients OA2qSLSM and OA2qSRSM. In a second step, these operators are re-normalized to the standard definition of the Wilson coefficients C 7 and C 7 . The initial coefficients OA2qSR, OA2qSL are not discarded, but co-exist besides CC7, CC7p. Thus, the user can choose in the implementation of the observables which operators are more suitable for his purposes.

Validation
The validation of the FlavorKit results happened in three steps: 1. Agreement with SM results: we checked that the SM prediction for the observables agree with the results given in literature 2. Independence of scale in loop function: the loop integrals for two and three point functions (B i , C i ) depend on the renormalization scale Q. However, this dependence has to drop out for a given process at leading order. We checked numerically that the sum of all diagrams is indeed independent of the choice of Q. 3. Comparison with other tools: as we have pointed out in the introduction, there are several public tools which calculate flavor observables mostly in the context of the MSSM. We did a detailed comparison with these tools using the SPheno code produced by SARAH for the MSSM. Some results are presented in the following.
We have compared the FlavorKit results using SARAH4.2.0 and SPheno3.3.0 with superiso 3.3 -SUSY_Flavor 1 and 2.1 -MicrOmegas 3.6.7 -SPheno 3.3.0 a SPheno version produced by SARAH 4.1.0 without the FlavorKit functionality Since these codes often use different values for the hadronic parameters and calculate the flavor observables at different loop levels, we are not going to compare the absolute numbers obtained by these tools. Instead, we compare the results normalized to the SM prediction of each code and thus define, for an observable X, the ratio X SM is obtained by taking the value of X calculated by each code in the limit of a very heavy SUSY spectrum.
As test case we have used the CMSSM. The dependence of a set of flavor observables as function of m 0 is shown in Fig. 7 and as function of M 1/2 in Fig. 8. We see that all codes show in general the same dependence. However, it is also obvious that the lines are not on as a function of m 0 using the FlavorKit (red), superiso (purple), SUSY_Flavor 1 (brown), SUSY_Flavor 2 (green), SPheno (blue), MicrOmegas (orange) and the old implementation in SARAH (red dashed). The three lines for SUSY_Flavor 2 correspond to different options of the chiral resummation. We used M 1/2 = 200 GeV, A 0 = 0, tan β = 10, µ > 0.
top of each other but differences are present. These differences are based on the treatment of the resummation of the bottom Yukawa couplings, the different order at which SM and SUSY contributions are implemented,  the different handling of the Weinberg angle, and the different level at which the RGE running is taken into account by the tools. Even if a detailed discussion of the differences of all codes might be very interesting it is, of course, far beyond the scope of this paper and would require a combined effort. The important point is that the results of FlavorKit agree with the codes specialized for the MSSM to the same level as those codes agree among each other. Since the FlavorKit results for all observables are based on the same generic routines it might be even more trustworthy than human implementations of the lengthy expressions needed to calculate these observables because it is less error prone. Of course, known 2-loop corrections for the MSSM which are implemented in other tools are missing. Finally, it is well known that the process B 0 s,d → ¯ has a strong dependence on the value of tan β. We show in Fig. 9 that this is reproduced by all codes.

Conclusion
We have presented FlavorKit, a new setup for the calculation of flavor observables for a wide range of BSM models. Generic expressions for the Wilson coefficients are derived with PreSARAH, a Mathematica package that makes use of FeynArts and FormCalc. The output of PreSARAH is then passed to SARAH, which generates the Fortran code that allows to calculate numerically the values of these Wilson coefficients with SPheno. The observables are derived by providing the corresponding pieces of Fortran code to SARAH, which incorporates them into the SPheno output. We made use of this code chain to fully implement a large set of important flavor observables in SARAH and SPheno. In fact, due the simplicity of this kit, the user can easily extend the list with his own observables and operators. In conclusion, FlavorKit allows the user to easily obtain analytical and numerical results for flavor observables in the BSM model of his choice.

A: Lagrangian
In this section we present our notation and conventions for the operators (and their corresponding Wilson coefficients) implemented in PreSARAH. Although a more complete list of flavor violating operators can be built, we will concentrate on those implemented in PreSARAH. If necessary, the user can extend it by adding his/her own operators.
The interaction Lagrangian relevant for flavor violating processes can be written as The first piece contains the operators that can trigger lepton flavor violation whereas the second piece contains the operators responsible for quark flavor violation. The general Lagrangian relevant for lepton flavor violation can be written as The first term contains the − − γ interaction, given by Here e is the electric charge, q the photon momentum, P L,R = 1 2 (1 ∓ γ 5 ) are the usual chirality projectors and α,β denote the lepton flavors. For practical reasons, we will always consider the photonic contributions independently, and we will not include them in other vector operators. On the contrary, the Z-and Higgs boson contributions will be included whenever possible. Therefore, the − − Z and − − h interaction Lagrangians will only be used for observables involving real Z-and Higgs bosons. These two Lagrangians can be written as where p is the β 4-momentum, and The general 4 4-fermion interaction Lagrangian can be written as where α,β,γ,δ denote the lepton flavors and Γ S = 1, Γ V = γµ and Γ T = σµν . We omit flavor indices in the Wilson coefficients for the sake of clarity. This Lagrangian contains the most general form compatible with Lorentz invariance. The Wilson coefficients A S LR and A S RL were included in [57], but absent in [42,58]. As previously stated, the coefficients in Eq.(A.6) do not include photonic contributions, but they include Z-boson and scalar ones. Finally, the general 2 2q four fermion interaction Lagrangian at the quark level is given by Here dγ denotes the d-quark flavor.
Let us now consider the Lagrangian relevant for quark flavor violation. This can be written as (A.10) The first two terms correspond to operators that couple quark bilinears to massless gauge bosons. These are Here T a are SU (3) matrices. The Wilson coefficients Q L,R 1,2 can be easily related to the usual C ( ) 7,8 coefficients, sometimes normalized with an additional 1 16π 2 factor. The 4d four fermion interaction Lagrangian can be written as where d α,β,γ,δ denote the lepton flavors. Again, we omit flavor indices in the Wilson coefficients for the sake of clarity. The 2d2 four fermion interaction Lagrangian is given by Here γ denotes the lepton flavor. L 2d2 should not be confused with L 2 2d . In the former case one has QFV operators, whereas in the latter one has LFV operators. This distinction has been made for practical reasons.
The 2d2ν and du ν terms of the QFV Lagrangian are Note that we have not introduced scalar or tensor 2d2ν operators, nor tensor du ν ones, and that lepton flavor (denoted by the index γ) is conserved in these operators. Finally, we have also included a term in the Lagrangian accounting for operators of the type (dΓ d)S and (dΓ d)P , where S (P ) is a virtual 10 scalar (pseudoscalar) state. This piece can be written as

B: Operators available by default in the SPheno output of SARAH
The operators presented in Appendix A have been implemented by using the results of PreSARAH in SARAH.
Those are exported to SPheno. We give in the following the list of all internal names for these operators, which can be used in the calculation of new flavor observables.

B.1: 2-Fermion-1-Boson operators
These operators are arrays with either two or three elements. While operators involving vector bosons have always dimension 3 × 3, those with scalars have dimension 3 × 3 × ng. ng is the number of generations of the considered scalar and for ng = 1 the last index is dropped.
These operators are derived by PreSARAH with the following input files    18 19 O u t p u t F i l e = "Gamma2Q .m" ;    16 17 O u t p u t F i l e = "Gluon2Q .m" ; 18 19 F i l t e r s = { } ; The normalization is changed to match the standard definitions by These operators are derived by PreSARAH with the following input files     18 19 O u t p u t F i l e = "Gamma2l .m" ; 20 21 F i l t e r s = { } ; The normalization is changed to match the standard definitions by In the following we omit flavor indices for the sake of simplicity. These operators are derived by PreSARAH with the following input files    16 17 O u t p u t F i l e = " Z 2 l .m" ; 18 19 F i l t e r s = { } ; O u t p u t F i l e = " H2l .m" ; 17 18 F i l t e r s = { } ; (dΓ d)S and (dΓ d)P

Variable
Operator Name Variable Operator Name These auxiliary 11 operators are derived by PreSARAH with the following input files   19 20 O u t p u t F i l e = "H2q .m" ;  19 20 O u t p u t F i l e = "A2q .m" ; We will denote the 4-fermion operators involving two leptons and two down-type quarks depending on whether they lead to LFV or to QFV processes: dd for LFV and dd for QFV.

Here we have definedÂ
(C.20) The mass of the leptons in the final state has been neglected in this formula, with the exception of the dipole terms K L,R

2
, where an infrared divergence would otherwise occur due to the massless photon propagator. Eq.(C. 19) is in agreement with [58], but also includes the coefficients A S LR and A S RL .
We give also here a description of the implementation of the different observables using the operators present in the SPheno output of SARAH.

97
Real ( dp ) , I n t e n t ( i n ) : : xc , x t C.2.9: P → ν Although P → ν, where P = qq is a pseudoscalar meson, does not violate quark flavor, we have included it in the list of observables for practical reasons, as it can be computed with the same ingredients as the QFV observables. The decay width for the process P → αν is given by [139] Γ (P → αν) = |G F f P (m 2 Here f P is the meson decay constant, mq and m q are the masses of the quarks in the meson and the Wilson coefficients G I XY are defined in Eq.(A. 16). The sum in Eq.(C.102) is over the three neutrinos (whose masses are neglected).