FlexibleSUSY extended to automatically compute physical quantities in any Beyond the Standard Model theory: Charged Lepton Flavor Violation processes, Higgs decays, and user-deﬁned observables

FlexibleSUSY is a framework for the automated computation of physical quantities (observables) in models beyond the Standard Model (BSM). This paper describes an extension of FlexibleSUSY which allows to deﬁne and add new observables that can be enabled and computed in applicable user-deﬁned BSM models. The extension has already been used to include Charged Lepton Flavor Violation (CLFV) observables, but further observables can now be added straightforwardly. The paper is split into two parts. The ﬁrst part is non-technical and describes from the user’s perspective how to enable the calculation of predeﬁned observables, in particular CLFV observables. The second part of the paper explains how to deﬁne new observables such that their automatic computation in any applicable BSM model becomes possible. A key ingredient is the new NPointFunctions extension which allows to use tree-level and loop calculations in the model-independent setup of observables. Three examples of increasing complexity are fully worked out. This illustrates the features and provides code snippets that may be used as a starting point for implementation of further observables.


Introduction
Exploring the parameter space of Beyond the Standard Model (BSM) theories, researchers frequently employ software packages to automate both intricate calculations and parameter scans.Nevertheless, these software tools are often restricted to a limited range of models, closely related to the Standard Model (SM), the Minimal Supersymmetric Standard Model (MSSM), the Two-Higgs-Doublet Model (2HDM), and their extensions involving higher-dimensional operators.This is only a slice of all intriguing possibilities that include a broad set of models.
To address this limitation, FlexibleSUSY [1,2] was created to be used for a broad class of supersymmetric (SUSY) or non-SUSY models, providing a tool for the comprehensive investigation of diverse theoretical scenarios. 1It is a software application primarily implemented in the Wolfram Language [7] and C++ based on SARAH [3,[8][9][10] and components from SoftSUSY [11,12], designed to produce an efficient and accurate C++ spectrum generator (a program searching for a consistent set of model parameters and calculating the pole mass spectrum and a set of observables) for physical models specified by the user.
The produced C++ program applies user-defined boundary conditions at up to three distinct energy scales within the model, incorporating Renormalization Group Equation (RGE) evolution between these scales.It further generates a collection of mixing matrices, pole masses, and auxiliary quantities.Recent versions of FlexibleSUSY have introduced the computation of several important observables that are suitable for phenomenological investigations and comparisons with experimental data.In particular, we highlight the extensions FlexibleAMU and FlexibleCPV introduced in Ref. [2] (they are responsible for the calculations of the anomalous magnetic moment and Electric Dipole Moment (EDM) of leptons), an update for FlexibleMW on precise calculation of the W -boson pole mass from Ref. [13], and FlexibleDecay [14] (a tool to calculate decays of scalars in a broad class of BSM models).
A key point of these observables is that they are integrated in FlexibleSUSY such that they are ready to be computed for any desired BSM scenario FlexibleSUSY is applied to.In order to achieve this, the mentioned extensions store information about the observables in suitable Wolfram Language and C++ meta code which is then automatically converted into actual C++ code specifically for each BSM scenario.
So far, new observables were added by individually modifying the internals of FlexibleSUSY, hence users could not add new observables in such a model-independent way.In the present paper, we explain a new FlexibleSUSY 2.8 design structure which solves this problem.It allows to integrate new observables on the meta code level, and it provides powerful options to define and finetune the computation of new observables without having to touch internals of FlexibleSUSY.
A number of new observables has already been integrated by using this new structure (they correspond to various Charged Lepton Flavor Violation (CLFV) processes), and in the future, further additional observables may be integrated either by FlexibleSUSY developers or by users of FlexibleSUSY.
To streamline and modularize the integration of new observables into FlexibleSUSY, an extension named NPointFunctions [15] was developed.This extension serves to automate the calculation of amplitudes and other quantities that rely on them for any high-energy model supported by FlexibleSUSY.In essence, NPointFunctions incorporates a well-defined approach of widely used packages FeynArts [16], FormCalc [17], and ColorMath [18] into FlexibleSUSY (up to technical implementation details to be mentioned later in appropriate sections).
The article is separated into two parts depending on the readers' interests: 1. Section 2 describes how to install and use FlexibleSUSY to calculate any of the available observables.This section is of interest for all users of FlexibleSUSY who may want to switch on the computation of observables.Reading it does not require knowledge of the internal structure of FlexibleSUSY.To get physical insights about the new CLFV observables, one is invited to read Appendix A.
2. Section 3 presents details on how to implement new observables.It provides a general outline and background information relevant for all observables, and it covers three specific examples of increasing complexity.It thus illustrates the range of possibilities and equips users with code snippets which can be used as a basis for further developments.Interesting features improving the functionality of FlexibleSUSY are mentioned separately in Appendix B.

Available observables and how to calculate them with FlexibleSUSY
Observables that are currently available in FlexibleSUSY are shown in Table 1.This section shows the reader how to calculate them with FlexibleSUSY.

Installation
There are no changes to the previous FlexibleSUSY version of Refs.[2,14] with respect to mandatory installation steps.All missing dependencies will be highlighted by FlexibleSUSY during the execution of the configure script (their list and hints about the installation can be found at developer's repository, see the program summary).To use observables that rely on NPointFunctions module (in particular, CLFV ones), FeynArts and FormCalc must be installed.

Output defined observable with FlexibleSUSY
To switch on the calculation of a desired, predefined observable O i (for example, bsgamma or LToLConversion from Table 1) for a physical model M a (like SM or MSSM), one needs to modify Observable FlexibleSUSY name

Loop level
Hints and comments ∆a ℓ

EDM
Electric dipole moment of a lepton [2] See Table 2 and Appendix A.2 µ-e conversion LToLConversion 0-1 See Table 3 and Appendix A.3 b → sγ bsgamma 1 --FlexibleDecay known SM (0-4), LO (0-1) for BSM Decays of scalars [14]  either model _ files/M a /FlexibleSUSY.m or models/M a /FlexibleSUSY.m.These files define the C++ spectrum generator output in the SUSY Les Houches Accord (SLHA) [19,20] or the Flavour Les Houches Accord (FLHA) [21] formats.The first file mentioned above is used by the script createmodel (see Sections 3.7.1-3.7.3) to create both the directory models/M a / and the second file, while the latter is executed by commands configure and make that create C++ spectrum generator itself.This means, that to always include the desired observables after the directory models/M a / is purged or cleaned, one modifies the first file.
The required modification of the file FlexibleSUSY.m is in the list ExtraSLHAOutputBlocks.Each observable that should be computed and appear in the spectrum generator output needs a corresponding entry which specifies the details of the observable and how it should appear in the output.For example, µ-e conversion is switched on in the Minimal R-Symmetric Supersymmetric Standard Model [22] (MRSSM) (for its particular FlexibleSUSY configuration called  The numerical value of µ-e conversion after the execution of the spectrum generator will be stored in SLHA format under the user-chosen number 41 in the block FlexibleSUSYLowEnergy (see meta/WriteOut.m).More details for arguments of LToLConversion are provided in Table 3.
To numerically calculate all added observables (apart from FlexibleDecay), the observable calculation must be enabled in the SLHA input file for C++ spectrum generator: In templates/observables/ @O i _ filename@.hpp.in@O i _ filename@.cpp.in In models/M a /observables/ M a _ @O i _ filename@.hppM a _ @O i _ filename@.cppThe coefficients are defined in files meta/Observables/O i /WriteOut.m.

File structure of new observables O i
Two new FlexibleSUSY features are described in this paper: a way to add new observables to FlexibleSUSY (implying new auxiliary files), and a way to generate C++ code to numerically calculate amplitudes that are required by these observables (NPointFunctions extension).This section addresses both extensions and explains how to implement new observables with the help of toy examples and code snippets of already implemented observables.
To create a spectrum generator, FlexibleSUSY usually first runs Wolfram Language code (located in meta/) to obtain model-specific expressions for mass matrices, self-energies, amplitudes, etc.The obtained expressions are converted to C++ form and are filled into C++ template files (located in templates/) to generate model-specific C++ code (located in model/M a /). Figure 1 shows the files relevant for the calculation of observables.The top left block shows Wolfram Language files in meta/ that contain general routines and observables implemented in previous versions of FlexibleSUSY.The top right block shows the corresponding observable-specific Wolfram Language files (located in meta/Observables/O i /) and the necessary C++ template files (located in templates/observables/).The bottom block shows the generated relevant C++ files that are eventually compiled and combined with other C++ files to the final spectrum generator executable (models/M a /run _ M a .x).
In the following, by creating a new observable we mean the generation of a specific C++ function that calculates the observable within the final spectrum generator.That is, one needs to create several C++ code blocks in appropriate places of the generated C++ spectrum generator, as the following template-listing with observable-and model-specific blocks shows: 3Listing 4: C++ template code to calculate an observable for a specific model Ma.
// 1. Function definition in Ma _ @Oi _ filename@.cpp:namespace Ma _ @Oi _ namespace@ { @Oi _ output _ type@ @Oi _ calculate(...)@ { // Complicated function body generated during the meta phase and specific for Oi and Ma } } // 2. Internal calculation of the observable: observables.@Oi_ name@ = Ma _ @Oi _ namespace@::@Oi _ calculate(...)@; // 3. Writing the observable in Les Houches format to an output stream: block << FORMAT _ * (observables.@Oi_ name@, ...); In the C++ source code snippet above, @O i _ name@ is replaced by the name of the observable and @O i _ output _ type@ is its C++ type.The @O i _ calculate(...)@ pattern will be replaced by the name of function that performs the calculation, which is defined in the M a _@O i _ namespace@ namespace.The body of the function will depend on the observable and the model under consideration.
In the following subsections it is described how the expressions are determined that replace the different patterns in the above C++ template code snippet.The files where these expressions are determined are the Wolfram Language files located in the directory meta/Observables/O i /, which will be described in the following subsections.In general, to define a new observable one creates five new files filled with the content described in the following steps: 1.In Section 3.1 we show how to define the C++ name of the observable, the observable's C++ type, etc. (done in the file meta/Observables/O i /Observables.m).
2. In Section 3.2 we show how to connect the calculation of the observable to the spectrum generator output (done in the file meta/Observables/O i /WriteOut.m).
3. In Sections 3.3-3.4we describe how to fill the body of the function that calculates the numeric value of the observable for a given parameter point (done in the meta-phase file meta/Observables/O i /FlexibleSUSY.m in Section 3.3 and two C++ template files in Section 3.4).
4. As a last step, we show how to modify the model-specific FlexibleSUSY model file named models/M a /FlexibleSUSY.m to register the observable to the model under consideration.
Each subsection will begin with general explanations and definitions.Then three recurring examples will be used for illustration.The first example illustrates the definition of an observable with minimal complexity -the output of a single numerical constant, where the numerical constant is not hard-coded but can be specified separately for each model FlexibleSUSY is applied to.The second example outputs the value of fermion masses, where the decision of which fermions to select can be specified separately for each model.The third example outputs the value of one-loop self-energies and thus illustrates how to use loop calculations in the definition of observables.

General case
The creation of a new observable O i starts with the file meta/Observables/O i /Observables.m.This file is supposed to define the C++ tokens such as the observable name (@O i _ name@) or the observable's C++ type (@O i _ output _ type@), etc.The file contains only one call to the function DefineObservables, which defines all C++ tokens.The general structure of this function call is shown in the following Wolfram Language source code snippet: Listing 5: General content of meta/Observables/Oi /Observables.m.
All further arguments of the function DefineObservables define how C++ tokens (for example, @O i _ output _ type@) for the observable O i should be replaced in code from Listing 4. Note that the names of the parameters defined for the observable (like parA, parB, etc.) are replaced in the strings by their values given in the model-specific FlexibleSUSY model file (the only exception: arguments in the prototype).For example, if one writes in the FlexibleSUSY model file FlexibleSUSYObservable O i [Fe,3], then the strings parA and parB are replaced by Fe and 3, respectively, in all strings on the r.h.s. in the function DefineObservable, so that for example the execution of the function GetObservableName produces "name _ with _ Fe3". 4  The meaning of two other mandatory arguments of the function DefineObservables is the following: GetObservableType defines the C++ type of the observable and replaces @O i _ output _ type@.In general, it can be any (real or complex) scalar, array or matrix C++ type.The syntax {1} in Listing 5 defines a complex-valued array of length 1 (corresponds to the C++ type Eigen::Array<std::complex<double>,1,1>).See also BrLToLGamma/Observables.m for another example.Note that array, vector or matrix types are convenient for storing Wilson coefficients, debugging information, etc.As described later, the connection between the observable type and the FlexibleSUSY Les Houches output is specified in the file named meta/Observables/O i /WriteOut.m.
GetObservablePrototype defines the C++ prototype for the function that calculates the observable (@O i _ calculate(...)@).Note that the function name is modified with the values of the parameters of the first argument of the function DefineObservable (parA is replaced with value of parA _ in Listing 5), while the function arguments are kept intact (parB is kept). 5  Let us now illustrate how to create new observables in FlexibleSUSY starting from the usage of the function DefineObservable in Observables.mfile.The examples are of increasing complexity and they will be continued in the next subsections.

Example 1: output a given number
Let us create an observable, called ExampleConstantObservable, that outputs a numerical constant specified by the user in the FlexibleSUSY.mfile, so that the value of the numerical constant is not hard-coded in the body of the observable C++ code but rather configured with the Wolfram Language files.This example will demonstrate the basic usage of the function that calculates an observable and allows us to familiarize ourselves with the workflow, as no physical calculation is required.First, one creates the meta/Observables/ExampleConstantObservable/ directory.Afterwards one creates the file Observables.minside this directory with the following content: 4 More advanced ways to name the C++ code parts are supported: 1.One may use expression like "$(parA+1)", which evaluates the content inside the parenthesis after substitution of observable pattern values (parA _ ). 2. The description of observable in the SLHA output can be generated from its name automatically by replacing underscores with spaces.This can be redefined by GetObservableDescription. 3.Both names for C++ templates and C++ namespace can be automatically generated from the Wolfram Language name of the observable by inserting an underscore before capital letters (or before numbers in front of them) and lowering the letter case.If the generated replacement for @Oi _ namespace@ (@Oi _ filename@) leads to undesired side effects, then one may use GetObservableNamespace (GetObservableFileName) to override the default behavior.4. GetObservableName uniquely defines the C++ name of the observable, which replaces @Oi _ name@.This name is used to store the numerical value of the observable internally and is generated automatically but can be overridden. 5The following two additional function arguments may be added: auto model and auto qedqcd.The first one allows to access all numeric quantities of the model (model parameters, masses, etc.).The second one provides access to known low-energy input observables.As in the general case discussed above, the function call of DefineObservables defines the Wolfram Language symbol of the observable to be ExampleConstantObservable, and it specifies that the observable depends on one parameter num _ ; the meaning of the parameter will be the value of the numerical constant to be output.Accordingly, as C++ return type for the function that calculates the observable we use an array of size one in GetObservableType.GetObservablePrototype defines the prototype of the function that performs the calculation and we use double num as function argument.

Example 2: show fermion masses
Let us turn to a little more involved example to illustrate the usage of model-specific quantities in the calculation of an observable and the possibility of filling the C++ template files with modelspecific content.As an example, we define an "observable" which can output several masses whose values are determined in FlexibleSUSY.Specifically we choose to allow the output the Modified Minimal Subtraction Scheme (MS)/Dimensional Reduction Scheme (DR) mass of a user-selected fermion at a given renormalization scale and/or the pole mass of a SM lepton.Again, similar to the example above, the fermion field and its generation number will be specified by the user in the FlexibleSUSY.mfile, so that these values are not hard-coded in the body of the observable C++ code but are configured with the Wolfram Language files.
To set up the observable, we create the directory meta/Observables/ExampleFermionMass/ and the file Observables.mwithin it, with the following content: Again, the function call of DefineObservables first defines the Wolfram Language symbol of the observable to be ExampleFermionMass, and it specifies that the observable depends on two arguments, merged into one Wolfram Language function.We assume that the concrete name of the fermion to output will be defined in the FlexibleSUSY model file models/M a /FlexibleSUSY.m, and we assume that it has the form fermion[gen], where fermion is the name of the fermion multiplet and gen is its index in the multiplet.Therefore, the observable ExampleFermionMass has the form of a function with the fermion _ [gen _ ] pattern as argument, which matches the multiplet name with fermion _ and the index in the multiplet with gen _ . 6We would like to return two numerical values, corresponding to the MS/DR mass of the fermion fermion[gen] at a given renormalization scale and the pole mass of a SM lepton of generation gen.Hence we define the observable type to be an array of length 2 in GetObservableType.The prototype of the function that calculates the observable is defined in GetObservablePrototype.The function takes the index of the fermion in the multiplet as first argument.The remaining two arguments model and qedqcd allow the access to the model parameters, masses, mixing matrices, etc.In the name of the function prototype ex _ fermion _ mass, the string fermion is replaced by the value of the variable fermion, which may be Fe, Fd, etc. depending on the user-selected model M a .As we plan to calculate different masses depending on the SLHA output block (see Section 3.2.2),let us implement the description of the observable as specified with GetObservableDescription above to have an explicit indication of the generated numerical values.There, the strings fermion and gen are replaced by the user-specified values of fermion and gen, respectively.

Example 3: lepton self-energy with NPointFunctions
As a third example, let us implement an observable to calculate the self-energy of a userspecified lepton at the one-loop level.This example will illustrate how to use the NPointFunctions extension in both simple and advanced ways, as well as how to output quantities in the FLHA format.As in the previous examples, we start by defining the observable and a suitable C++ function prototype that allow us to select the lepton (from a multiplet) whose self-energy shall be calculated and the contributions that should be taken into account: We call the observable ExampleLeptonSE, which takes the multiplet name (field _ ) and the index of the particle in the multiplet (gen _ ) as first argument, as in the previous example.Furthermore, we'd like to be able to select a subset of the full one-loop contribution.To achieve this, we provide a second parameter (contr _ ) which we will use later to only calculate a certain part of the full one-loop self-energy.In general, a fermion self-energy can be split into different covariants involving left-handed or right-handed projection operators P L,R and possibly the covariant / p, where p µ is the external fermion momentum.Let us in this example for simplicity output only the coefficients of the two covariants / pP L and / pP R , so we choose an array of length 2. The C++ function needs the index of the particle in the multiplet (gen) and the model object (model) as arguments.

Content of the file O i /WriteOut.m General case
After the definition of the observable name, type, etc. in the file O i /Observables.m,one can now connect the numerical value(s) of the observable to the Les Houches output of the C++ spectrum generator in the file meta/Observables/O i /WriteOut.m.Currently, there are two automated ways to output the results of calculations in FlexibleSUSY: via the SLHA (exists since FlexibleSUSY 1.0) and the FLHA (the usage is explained in this article) formats.Both are defined in the meta/WriteOut.m file, which provides functions to write numbers, arrays or matrices to specified output blocks.
In the simplest case the value of an observable O i [parA, parB, ...] is a single complex number (see Examples 1-2 below or the ℓ i → ℓ j γ decay; their Les Houches output is shown in Sections 3.7.1-3.7.2).If we want to write its real part to the FlexibleSUSYLowEnergy SLHA output block, we define a function called WriteObservable in meta/Observables/O i /WriteOut.mas follows: The first argument to the WriteObservable function is a string with the SLHA output block name (here: "FlexibleSUSYLowEnergy"), which reflects the SLHA block name in the FlexibleSUSY model file models/M a /FlexibleSUSY.m.Since the observable type is a complex-valued array of length 1 in this example, we have to output the 0-th element of the array (in C++ index convention).We obtain the real part of the array element by applying the Re function.
As another example we consider the output of the real parts of several Wilson coefficients (see also Example 3 below and its Les Houches output in Section 3.7.3).This requires a more involved definition of the function WriteObservable: { "fermions1, operator1, Oalpha1, Oalphas1, contributions1, num _ value, \"comment1\"", "fermions2, operator2, Oalpha2, Oalphas2, contributions2, num _ value, \"comment2\"", ... }, { "num _ value" -> "Re(observables."<> Observables GetObservableName[obs] <> ")", ... } ]; We define the function WriteObservable to write the real parts of the Wilson coefficients to the FWCOEF Les Houches output block. 7The return value of WriteObservable is a list of strings, with each string consisting of a comma-separated tuple of a fermion name (fermions 1 ), an operator name (operator 1 ), etc., which should be replaced by appropriate values.See Ref. [21] for a description of the FLHA format.The numbering convention for the FLHA format in FlexibleSUSY is the following: Wilson coefficients must occupy the positions from the end of the C++ output array.In the definition of the function WriteObservable the substring "num _ value" is replaced by the real parts of the Wilson coefficients (accessed via GetObservableName) via the function StringReplace.Note that there is no need to explicitly specify the numbering of the Wilson coefficients in the replacement rule for "num _ value", as it is done automatically.

Example 1: output a given number
We continue our minimal Example 1 from Section 3.1.1,where the observable is defined to just be a user-defined numeric constant specified in the FlexibleSUSY model file.The second step is to connect the numeric value of the observable to the SLHA (or FLHA) format, which is the output of the spectrum generator.This is done by the following definition placed in the file As the first argument of the function WriteObservable shows, the numeric value of the observable is written to the SLHA output block FlexibleSUSYLowEnergy, which matches a corresponding definition in the model file models/M a /FlexibleSUSY.m.The function returns a string (which must be valid C++ syntax), where we take the real part of the first entry from the array of length 1, where the observable is stored (with the zero-based index convention in C++).Since for the output no information about the observable is needed, one uses the _ pattern in the specification of ExampleConstantObservable. The code listing above leads to the usage of the C++ parser for the SLHA format named FORMAT _ ELEMENT, see Listing 4.

Example 2: show fermion masses
Now we turn back to Example 2 from above and connect our observable (two given fermion masses) to the output of the spectrum generator.Similarly to Example 1, one creates the file meta/Observables/ExampleFermionMass/WriteOut.m and places the definitions of the function WriteObservable there.In this example we would like to output two different masses at the same time.To store them internally, we have defined the observable type to be an array of length 2, see the specification of GetObservableType in line 3 of Listing 7. To write the two fermion masses to the output, one could do the following: Here, we define two distinct behaviours of our observable with two different Les Houches output blocks (which must be reflected by appropriately named lists in models/M a /FlexibleSUSY.m): "FlexibleSUSYLowEnergy" and "ExampleLeptonMass".In both cases, definitions lead to the C++ parser named FORMAT _ ELEMENT, see Listing 4. The fermion mass in the first entry of our twocomponent observable array is written to the block FlexibleSUSYLowEnergy, while the second entry is written to the block ExampleLeptonMass.Note that the name of this second block is not a part of the official SLHA standard, but is a non-standard addition we use for this example.

Example 3: lepton self-energy with NPointFunctions
In Example 3, we aim to output the two components of a lepton self-energy.We can use the automatic FLHA output format to achieve this: In the definition of WriteObservable we provided a pattern for the index of the lepton in the lepton multiplet (gen _ ) explicitly in the first argument of the observable name ExampleLeptonSE, because the concrete value of this index must be known to select the appropriate self-energy for the output.

General case
In the previous sections it was described how to define a new observable (done in the file Observables.m)and how to write the numeric value of the observable to the Les Houches output of the spectrum generator (defined in the file WriteOut.m).In this section we describe how to generate the content of the C++ function, that calculates the numeric value of the observable.This C++ function is defined in the C++ template files located in templates/observables/, which contains placeholders that are replaced by expressions to calculate the observable.The rules that specify how to replace the placeholders by appropriate expressions are defined in the observable-specific file meta/Observables/O i /FlexibleSUSY.m.This file contains at least the function WriteClass, which performs the following tasks: 1. Fill the C++ templates from templates/observables/ with appropriate C++ code.
2. Move the filled C++ templates into models/M a /observables/.

} ];
The function WriteClass has three parameters: the explicit name of the observable (obs), the list of all observables that are requested -by the variable ExtraSLHAOutputBlocks from the file models/M a /FlexibleSUSY.m -to be calculated (allObs _ ), and a list of the C++ template file names, where the replacements should be performed, and the corresponding output file names (files _ ).The body of the function WriteClass consists of three parts: 1.The first part is the If statement, where all C++ expressions are created (Task In the above example source code listing prototypes strings are created by replacing the "@type@" and "@prototype@" tokens by the observable's C++ type and the prototype of the function that calculates the numeric value of the observable, respectively. 2. In the second part (Task 2 indicated in the listing) the C++ tokens ("@npf _ headers@", "@npf _ definitions@", etc.) are replaced in the C++ template files (files) that are passed to the function with the help of the ReplaceInFiles function.
3. The third part (Task 3) is to specify the function's returned expression.In the example above the function returns a list of replacement rules, whose use is described below.
Note that in Task 2 the tokens in the C++ template files can be replaced by strings that can in principle be arbitrarily large.However, we recommend to put as much generic information as possible into the C++ template files and replace the tokens in the template files only by model-specific information.For the latter we recommend to use the full power of FlexibleSUSY's helper routines and functions located in meta/TextFormatting.m, meta/CConversion.m and meta/Utils.m.We refer the reader also to Section 3.5, where we describe the NPointFunctions extension and how one can use it to generate C++ code for amplitudes.
The expression returned by the WriteClass function (Task 3) should be a list of replacement rules, which gets stored internally in the variable ObservablesExtraOutput["O i "].The replacement rules stored in the returned list can be accessed and re-used later, if needed.For example, one can access expression stored in the rule defined in line 40 from Listing 14 as follows: Besides the possible manual re-use of the returned expressions, FlexibleSUSY automatically performs the following two tasks with the returned list of replacement rules: 1.If there exists a replacement rule of the form "C++ vertices" -> list 1 , then the expression list 1 must be a list of vertices that are required to calculate the observable.Each vertex is represented by a list of SARAH fields, e.g.{bar[Chi], Chi, VP} in the MRSSM (represents the χ0 a χ 0 b γ vertex with two neutralinos and a photon).For each required vertex FlexibleSUSY creates a corresponding C++ function for its numerical evaluation when calculating the numeric value of the observable.2. If there exists a replacement rule of the form "C++ replacements" -> list 2 , then the expression list 2 must be a list of replacement rules that should be applied to all C++ template files, see Appendix B.3 for an example.

Enabling optional NPointFunctions extension, simple settings
If the observable relies on tree-level or one-loop level Feynman diagrams then there is an automated way to generate them in FlexibleSUSY with the help of the NPointFunctions extension already announced in Ref. [15].NPointFunctions is typically used from within the WriteClass function.Internally, the extension calls FeynArts [16], FormCalc [17], and ColorMath [18] to generate analytic expressions and converts them into C++ form for their numeric evaluation in FlexibleSUSY.It thus allows access to Feynman diagrammatic computations in the definition of observables.
NPointFunctions can be used in two modes which we will refer to as simple and advanced modes.The modes differ by the accessible settings.Both types of settings serve to modify the calls of the FeynArts and FormCalc routines: In the simple mode only the settings listed in Table 4 are accessible, which allow topology-indepedent modifications.The advanced mode allows many more settings, enabling in particular to select specific option values for selected topologies.In the present section, we focus on the simple settings from Table 4.These simple settings are also illustrated with the help of Example 3 in Section 3.3.3.Later, the advanced settings are explained in detail in Section 3.5 and exemplified in an extra Section 3.5.7.
ZeroExternalMomenta can currently be True, False, OperatorsOnly, and ExceptLoops.This option also specifies the way external masses should be treated for specific topologies, fermions bispinors, and in loop integrals.If set to True, all external momenta are set to zero everywhere, while if set to ExceptLoops, the scalar products in loop integrals are kept.This allows one to correctly implement the expressions for self-energy-like diagrams relevant for CLFV processes in particular.
UseCache stores the FeynArts and NPointFunctions output for future usage, if enabled.This speeds up the code generation if model/M a / is purged or parts of the meta code are modified.
KeepProcesses specifies the mode of NPointFunctions.There are two modes to calculate an amplitude: using simple settings only (Observable -> None) and with the help of advanced ones defined in meta/Observables/O i /NPointFunctions.m files (Observable -> O i []), see Section 3.5. 8

Example 1: output a given number
Each observable undergoes various calculational steps before being evaluated into an actual numerical result.In general, the definitions of the calculational steps are distributed between the Wolfram Language file meta/Observables/O i /FlexibleSUSY.m and the C++ template files discussed.Specifically, the WriteClass function in the FlexibleSUSY.mfile is supposed to generate model-specific expressions to be placed into C++ template files from Section 3.4.1.As stated before, in general it is recommended to put as much information as possible into the C++ template files and keep the WriteClass function minimal.In case of Example 1, however, the example is so simple that we would like WriteClass to generate the complete C++ code for all calculations (just output a number in this case), thus the C++ template files will not contain much code.This code generation is done with the definition of @type@ @prototype@ { @type@ res {num}; return res; }", { "@type@" -> CConversion CreateCType[Observables GetObservableType[#]], "@prototype@" -> Observables GetObservablePrototype[#] } ] &/@ observables; ... ( * Task 3: return an empty list * ) {} 8 The simple mode with the simple settings discussed here was developed mainly to allow users a simple starting point.It is also used in unit testing of FlexibleSUSY.Currently, FlexibleSUSY performs unit tests of self-energy expressions in the SM/MSSM and Z-boson penguins in the MRSSM.It is also used in Example 3 to demonstrate the basic usage of the NPointFunctions package and can serve to write the first iterations of the code that relies on NPointFunctions.Typically, the code for actual physics observables such as the ones discussed in Section 2 will use the advanced settings of NPointFunctions.
In the above source code listing the variable definitions is defined to contain the entire definition of the function that calculates the observable (including function body).It is generated by replacing the @type@ token by the concrete C++ type Eigen::Array<std::complex<double>,1,1> and the function name @prototype@ (including its argument double num) as specified by the definitions of Section 3.1.1 in a generic string.In the body of the function, we initialize a local variable res of type @type@ and set its value to the value of num.The value of res is then returned to further fill the SLHA block entries.The final function definition in the definitions variable will be used later in Section 3.4.1 together with the C++ template files to construct the complete C++ code to output the numerical value of the observable.
The WriteClass function finally returns an empty list (Task 3), because we do not need any of the defined expressions anywhere else in this example.
Note that the function definition generated in Task 1 is later put into the C++ template files in templates/observables/ an is thus closely connected to their content, as shown in Section 3.4.1.

Example 2: show fermion masses
Now we return to Example 2 where we output fermion masses.In this example we follow the general recommendation to put as much C++ code as possible into the C++ template files and as little as possible into FlexibleSUSY.m.Thus, in this example we fill the following code into the WriteClass function in the FlexibleSUSY.mfile of Listing 14: @type@ @prototype@ { return forge<@type@, fields::@fermion@>(gen, model, qedqcd); }", { "@fermion@" -> SymbolName[Head[First[#]]], "@type@" -> CConversion CreateCType[Observables GetObservableType[#]], "@prototype@" -> Observables GetObservablePrototype[#] } ] &/@ observables; ... ( * Task 3: return an empty list * ) {} Again, the variable definitions contains the definition for each C++ function for the observable calculation of the function that calculates the numerical value of the observable.This function has the name and type that were specified via GetObservablePrototype and GetObservableType, respectively.In this example the function body contains only a single line, which contains a call of the template function forge.By this call we delegate the computation of the observable to the forge function, which is defined entirely at the C++ level in Section 3.4.2.The same strategy is applied in several of the predefined observables discussed in Section 2 such as ℓ i → ℓ j ℓ k ℓ c k or µ-e conversion.The template parameter fields::@fermion@ of the forge function above is an internal C++ type that represents the fermion whose mass shall be output (in the MRSSM, for example, the token @fermion@ is replaced by Fe).
Again, we do not want to use any expressions defined in the WriteClass function anywhere else, so we return an empty list from the function (Task 3).

Example 3: lepton self-energy with NPointFunctions (Observable -> None)
This example illustrates the use of the NPointFunctions extension of FlexibleSUSY.We will use the NPointFunctions extension in the simplified mode, where we do not make use of the advanced settings in NPointFunctions.m and just specify the option Observable -> None.The following code snippet shows the content of the WriteClass function: TextFormatting ReplaceCXXTokens[" @type@ @prototype@ { const auto npf = npointfunctions::@name@(model, {gen, gen}, {}); "@prototype@" -> Observables GetObservablePrototype[#], "@name@" -> name First of all, as NPointFunctions returns the same code for different generations of external particles, we remove potential duplicates by replacing integer generation with the _ pattern for simplicity (it could have been any symbol or number), see Task 1a.
Afterwards, we run NPointFunctions (Task 1b) to generate the relevant C++ code, and all required C++ vertices.For this we create a local function (defined as a Module), which we map to all observables in the observables list.Note, however, that in more involved calculations it may be better to define a non-local, separate function to process all observables, instead of a local one via a Module, as done here.The local function to obtain the C++ code to numerically evaluate the observable(s) stores its results in the following local variables: field represents the lepton particle that we specify in the definition of the observable (like Fe for example in the SM or MRSSM).
contr contains an expression that allows one to select certain contributions from the self-energy that should be output.In this example we do not make use of this selection feature.See Section 3.5.7 for an advanced example.
npf contains the main result of NPointFunctions: an object with generic amplitudes and, sometimes, replacement rules.See Table 4 for a list of all possible options.In this example, we calculate the Fe -> Fe amplitude, from which we will extract the self-energy (note that the outgoing particle is typed as Fe and not bar[Fe]).We do not store the results in the cache.OnShellFlag is used internally by the option FormCalc OnShell for the function FormCalc CalcFeynAmp.ZeroExternalMomenta is passed to the function FormCalc OffShell.The loop level is specified explicitly to be 1.The renormalization scheme is decided to be set by FlexibleSUSY.
We chose a simple operation mode of the NPointFunctions extension by setting Observable -> None.This implies, that no optional file O i /NPointFunctions.m with advanced settings is used, see Section 3.5 for more details and the usage example in Section 3.5.7.The option KeepProcesses allows selecting certain topologies, provided by the list of default values within the function GetExcludeTopologies in meta/NPointFunctions/Topologies.m.
To extend this list, we refer to the usage of FeynArts $ExcludeTopologies.
name contains the C++ name of NPointFunctions-generated function that calculates the numerical value of the observable.
basis is a list of replacement rules to extract certain sub-expressions of the amplitude and store them in local C++ variables.The rhs. of each replacement rule defines the expression, whose prefactor should be extracted.The lhs. defines the name of the C++ variable in which this prefactor shall be stored.Note that patterns are currently not supported to select the subexpression.The replacement rules stored in the basis variable should be passed to the InterfaceToMatching function, which performs the extraction of the prefactors.9 cxxVertices contains the list of all vertices that are required to numerically calculate the observable.They have to be returned from the WriteClass function (see Task 3), because FlexibleSUSY stores the C++ code for the vertices in other files.
npfDefinitions contains strings with all generated amplitude-related C++ functions required for the numerical evaluation of the amplitude.The function CreateCXXFunctions converts the NPointFunctions object into C++ code.The user provides the name of the C++ function and specifies which color structures are expected (Identity or SARAH Delta).
definitions contains the C++ function that calculates the numerical value of the observable with the help of C++ code created by the NPointFunctions extension.The function body consists of a call to the generated irr _ se function (the function name is stored in the name variable), to which the model parameters (model object) and the indices of the incoming and outgoing particles (both are gen here) are passed.The last argument contains the set of momenta of the external particles.In the current implementation external momenta are assumed to be zero and so the last function argument should be an empty initializer list {}.
Since basis contains two replacement rules that define prefactors of certain Lorentz structures, irr _ se will return a two-valued array.The generated function eventually returns the two components of this array.
In Section 3.5.7 the example described here is extended to use advanced NPointFunctions options.

General case
The function WriteClass defined in the file meta/FlexibleSUSY.m in the previous section fills C++ templates with model-specific information, and places the resulting files into the directory models/M a /.The content of the C++ template files is discussed in this section.
The numerical calculation of the observable might require more C++ auxiliary functions or classes.To use extra functions/classes one could add appropriate preprocessor #include directives.
Note again, that there is some freedom to move code between the O i /FlexibleSUSY.m and @O i _ filename@.cpp.infiles.However, as stated before, we recommend to move as much C++ code as possible into the dedicated directory templates/observables/, as C++ source code is often more robust compared to Wolfram Language scripts due to static type checking.

Example 1: output a given number
To output a single number we create two C++ template files in templates/observables/, as stated above.They will be filled by the WriteClass function defined in FlexibleSUSY.mcreated in Section 3.3.1.The content of the header template file example _ constant _ observable.hpp.in is the same as in general Listing 19.In principle, we can fill the C++ definitions template file with the default content provided in Listing 20.Nevertheless, as we do not use the NPointFunctions module for the observable calculation (realized with lines 3-6 in Listing 16), we can simplify the file as follows: Listing 21: Content of C++ definitions template example _ constant _ observable.cpp.in.#include "@ModelName@ _ @filename@.hpp"namespace flexiblesusy::@namespace@ { @calculate _ definitions@ } // namespace flexiblesusy::@namespace@ In this file, the template C++ token @calculate _ definitions@ will be replaced with the content of the variable definitions defined in the file O i /FlexibleSUSY.m described in Section 3.3.1.In this way, during the meta phase of FlexibleSUSY the complete C++ code will be generated for the function calculate _ example _ constant _ observable (defined in Listing 6) which outputs the number num.We do not need to provide any additional observable-specific code, as everything is already provided by the WriteClass function in O i /FlexibleSUSY.m during the Wolfram Language meta phase.

Example 2: show fermion masses
As stated above, we need to provide two C++ template files in templates/observables/.They will be filled during the FlexibleSUSY meta phase and embedded into the rest of the C++ spectrum generator.In this example the content of the template header file example _ fermion _ mass.hpp.in is the same as in the general Listing 19.The template file example _ fermion _ mass.cpp.in,however, contains two changes compared to the generic example from Listing 20.First of all, all commands related to NPointFunctions can be omitted, since we do not use Feynman diagrammatic calculations in this example (like in Example 1).Apart from this, our goal with this example is to demonstrate how one can move as many calculations as possible to the C++ template file, while keeping the flexibility to insert model-specific information.As discussed in Section 3.3.2,we prepared for this by having a minimal function body defined on the Wolfram Language level in the file ExampleFermionMass/FlexibleSUSY.m,where the function body only calls another C++ function named forge.This function shall now be defined.The following source code listing shows the content of the C++ template file that contains this definition: 10Listing 22: Content of C++ definitions template example _ fermion _ mass.cpp.in.#include "@ModelName@ _ mass _ eigenstates.hpp"#include "cxx _ qft/@ModelName@ _ qft.hpp" #include "@ModelName@ _ @filename@.hpp"#include "error.hpp"namespace flexiblesusy { using namespace @ModelName@ _ cxx _ diagrams; namespace @namespace@ { template <typename RTYPE, typename FIELD> auto forge(int idx, const @ModelName@ _ mass _ eigenstates& model, const softsusy::QedQcd& qedqcd) { context _ base context {model}; auto context _ mass = context.mass<FIELD>({idx});std::complex<double> lepton _ mass; switch (idx) { case 0: lepton _ mass = qedqcd.displayPoleMel();break; case 1: lepton _ mass = qedqcd.displayPoleMmuon();break; case 2: lepton _ mass = qedqcd.displayPoleMtau();break; default: throw OutOfBoundsError("fermion index out of bounds"); } RTYPE res {context _ mass, lepton _ mass}; return res; } @calculate _ definitions@ } // namespace @namespace@ } // namespace flexiblesusy In the code listing above, the template C++ token @calculate _ definitions@ will be replaced by the content of the variable definitions, defined in the ExampleFermionMass/FlexibleSUSY.m file in Section 3.3.2,which calls the forge function.The forge function defined above illustrates how to access two different types of particle masses.First, the running MS/DR mass of the fermion specified by FIELD and idx (which can correspond to, e.g.Fe [2] or Fd [3] at the Wolfram Language level in the selected BSM model) is obtained by calling the mass function template.Afterwards, the pole mass of the SM lepton of generation idx is obtained from the qedqcd object.The forge function finally returns an array of length 2 containing these two masses.

Example 3: lepton self-energy with NPointFunctions
As in the examples above, we have to create two C++ template files that are responsible for the numerical calculation of the observable.This example uses NPointFunctions to generate analytical expressions for the self energies, but it is otherwise a structurally simple example.Hence, like in Example 1 (but unlike Example 2) we do not delegate computations to a forge function; rather, the main definition of the observable is done by the function WriteClass in FlexibleSUSY.mcreated in Section 3.  We recall that in this example we use NPointFunctions in the simple mode specified by the option Observable -> None; the example will be continued in Section 3.7.3,which shows how to enable and use the computation of the self-energy in a concrete spectrum generator.An alternative version of the example is provided in Section 3.5.7,where the example is modified to illustrate the use of the advanced settings of NPointFunctions.

Content of the optional file O i /NPointFunctions.m, advanced settings
The NPointFunctions extension allows to generate Feynman diagrammatic calculations to obtain analytical expressions for observables in any specific model M a when the FlexibleSUSY meta phase is executed.The simple usage of NPointFunctions was demonstrated in Section 3.3.3,but frequently it is necessary to fine-tune calculations done with NPointFunctions.Such a fine-tuning is possible via advanced settings of NPointFunctions explained in this section.It may involve the selection of subsets of Feynman diagrams, remove contributions, selecting the regularization scheme, or specifying the fermion order in e.g.four-fermion amplitudes.Parts of the fine-tuning may be accessible to users of the observable (such as the arguments Vector, Scalars, etc. of observables in Tables 2-3), other parts may be found only in the definition of the advanced settings.
In the following we begin by explaining how to enable the advanced settings and giving an overview, then we will explain all advanced settings and finally provide an example.

Enabling advanced settings, overview
We begin by explaining how to enable and access the advanced settings when setting up the definition of an observable.The NPointFunctions extension is called from the observable-specific file O i /FlexibleSUSY.m discussed in Section 3.3.Listing 18 combined with Listing 14 provides an example for using NPointFunctions in its simple mode without advanced settings.To use the advanced mode, this file must be modified as follows: Here the line Observable -> O i [] enables the advanced settings (note that the variable obs was defined in the first line in Listing 14).This option requires the existence of the corresponding file meta/Observables/O i /NPointFunctions.m which should contain the detailed advanced settings described further below.
The second line defines the option KeepProcesses via the variable contr (instead of the hard-coded value Irreducible in Listing 18) and opens up an important way to access advanced settings.To explain the meaning of this line we briefly recall the overall structure of setting up an observable.From the user's perspective, an observable ultimately is called as described in Section 2, and observables may have options contr, see e.g.Tables 2-3, with possible values being a keyword or a list of keywords such as Vectors, Scalars, Boxes, etc.Using these options influences how the observable is evaluated.The definition of such keywords and how they influence the observable is part of the advanced settings of NPointFunctions.
Using such a contr option has to be prepared in the file O i /Observables.mas exemplified in Section 3.1.3by specifying the appropriate argument in the definition of the observable.Via line 9 of Listing 18 and line 3 of Listing 23 the variable contr propagates into the option of KeepProcesses.The construction in the code makes sure that the option is always a list.
The advanced settings are listed in Table 5 and explained in detail in the following Sections 3.5.2-3.5.6.As an overview, NPointFunctions relies on FeynArts and FormCalc to produce Wolfram Language expressions for model-specific class-level amplitudes, and the advanced settings fine-tune the calls to FeynArts and FormCalc routines, so familiarity with these tools is required to some extent.11Among the advanced settings, topologies, diagrams and amplitudes can be used to define the keywords which can be set by users as values of contr as explained above; these keywords then influence the execution of the functions FeynArts CreateTopologies, FeynArts InsertFields, and FormCalc CalcFeynAmp.The further advanced settings influence the calls to FeynArts and FormCalc and manipulate their output in ways which are fixed for the observable and can be influenced only by directly modifying the file O i /NPointFunctions.m.

topologies
The observable calculation by the NPointFunctions extension starts from the generation of topologies.It is possible to select particular topologies in the calculation by using the option KeepProcesses in line 3 of Listing 23 as explained above.The possible values of contr are (lists of) keywords: Vectors, Scalars, TreeLevelSChannel, etc.
In this section we focus on how to define such keywords, generally now called Contribution i , and to associate them with topologies of Feynman diagrams.(a) FeynArts topology names.The latter uniquely define topologies by connections to adjacency matrices, as shown in Figure 2.This explicit naming of topologies is useful, as they can also be used to modify diagrams and amplitudes, as described in the following subsections.
The name of a topology can be any Wolfram Language symbol that has to be connected with the actual adjacency matrix.All relevant topology names can be defined in a separate file meta/NPointFunctions/Topologies.m, as follows: Listing 24: Example for definitions of topologies (in meta/NPointFunctions/Topologies.m).

Special syntax
See the text description from Figure 2. One draws the FeynArts topology with explicit vertex numbers, creates an adjacency matrix, and eliminates all entries with redundant information: the zero topleft submatrix representing propagators between external particles, and entries below the diagonal as the adjacency matrix is symmetric.The rest is combined line by line into a one-dimensional list that is stored in meta/NPointFunctions/Topologies.m.These steps are done by the function AdjaceTopology defined there.
2. A connection between a set of simple topologies and their group name, like treeAll.This is useful to shorten the application of other settings.
If the topology names are defined like this, an example content of the file O i /NPointFunctions.m might define the keyword TreeLevelSChannel and associate it to the topology treeS as follows:

diagrams (generic level modifications)
The keywords Contribution i can also be associated with generic-level field patterns (in FeynArts nomenclature) to remove certain generic-level fields in the Feynman diagrams.Table 7 shows the required command in the file O i /NPointFunctions.m to define the associations with field patterns.The logic is the same as for topologies, but the main new ingredient is the Command b option.It defines the kind of generic fields that should be removed, and it should be a pure function which returns True or False and which takes several arguments while traversing the tree of FeynArts insertions: FeynArts TopologyList, the current topology, and the current class level insertion.Apart from standard Wolfram Language expressions, the following ingredients may be used to specify Command b : 12• Fields appearing in tree-parts of diagrams are specified by TreeFields.
• Fields appearing in loops are specified by LoopFields.
• Any generic field of scalar type FeynArts S, fermion type FeynArts F or vector boson type FeynArts V.
• A specific field, derived from external particles.For example, one can remove from the diagram fields that correspond to the external particles numbered 1 or 3 with the usage of the argument FieldPattern[#, 1|3] in the function FreeQ below.
Altogether, an example is given by the setting: Then, the function call NPointFunction[..., KeepProcesses -> {Vectors}] deletes the diagrams with the topology penguinT when the external particles 1, 3 do appear inside the loop.A second example, demonstrating the Absent mode of WhenToApply, is given by: In this way, if KeepProcesses -> {Vectors} is specified (thus Scalars is not specified), the third line of the listing applies and eliminates diagrams of penguinT topology which contain scalar particles in their tree-level parts.

amplitudes (class level modifications)
Often, modifications on the level of topologies and generic-level amplitudes allowed by the topologies and diagrams settings, together with other options of Table 4, provide sufficient finetuning.Sometimes, however, removing amplitudes on the classes-level is required.This is possible via specifying the amplitudes[LoopNumber, WhenToApply] setting.The syntax is identical to the one of diagrams and given in Table 7.As an example, amplitudes with massless particles can require special treatment, and the following setting allows to remove them from the generated expressions (assuming they will be handled properly elsewhere):

Abbreviation
Values Hints LoopNumber 0 or 1 As in Table 6 Momenta i ∀ Symbol from ZeroExternalMomenta See Table 4 Rule j

Special syntax See DiracChains.m
Table 8: Syntax and options for chains.Files in column "Hints" are located in meta/NPointFunctions/.

order, chains (Dirac algebra modifications)
NPointFunctions can be used to extract Wilson coefficients.For observables with external fermions, there might exist a need to change or fine-tune fermion chains to achieve the desired operator structure.This is done by the settings order and chains, where both the desired fermion order of external fermions as well as chains to be dropped are specified.
The syntax for order can be understood as follows.Imagine that we would like to obtain Wilson coefficients corresponding to the process µ − → e − e − e + , or more conveniently for the crossed 2 → 2 process µ − e − → e − e − instead, where the outgoing positron is replaced by an incoming electron.That means we need to obtain coefficients of expressions ū(e − )Γ i u(µ − ) ū(e − )Γ j u(e − ), where Γ represent structures involving Dirac matrices and momenta.For example, in the source code for the observable BrLTo3L, the NPointFunctions extension is accordingly called as where the variable lep will correspond to leptons during the meta phase execution for concrete models.Interpreting the incoming particles as muon and electron, we specify the required fermionic structure by the following line in the file NPointFunctions.m (the field numbers are as in FeynArts InsertFields): order[] = {3, 1, 4, 2}; Once the fermionic order is fixed by order, and amplitudes are calculated, one obtains multiple chains consisting of Dirac matrix products.Often certain Dirac chains may be simplified or neglected, thanks to the desired level of precision or the choice of basis of Wilson coefficients.Dropping specific types of Dirac chains can be implemented by using the chains setting described in Table 8.An example code which serves to neglect expressions proportional to the mass of the electron from the example above is given by: In general, the Rule j appearing on the right-hand side of this example (or the general case in Table 8) must be of the form ChainNumber[Entry1, Entry2, ...] -> 0 The syntax can be described as follows (the source code which interprets these expressions is in the file meta/NPointFunctions/DiracChains.m, and further details can be found there): ChainNumber is an integer which defines a chain number, e.g. for the fermion order {3,1,4,2} the chain numbered 1 consists of particles {3,1} and chain numbered 2 is {4,2}.
Entry 1 corresponds to a pattern for a Dirac chain which may be a non-commuting product of projection operators P L,R , γ µ matrices with open indices or γ µ matrices contracted with external momenta.We mimic FormCalc notation inside the Dirac chains so that the integer

Other modifications
Here we briefly describe and exemplify a set of settings which are typically of minor importance.
In some cases, one should not sum over all generations for some generic field in the amplitude.For example, in CLFV observables, there often appear penguin contributions with external selfenergy-like corrections.In such diagrams, the fermion in the tree-level propagator should differ from the external fermion.This behavior can be defined via the setting sum: This changes the expressions at the loop level LoopNumber = 1 in the following way.For the topology identified as inSelfT one modifies the summation over generic fields in the propagator under FeynArts number 6, so that the sum over the field being equal the first external particle is omitted on C++ level.
The setting momenta can be used to eliminate certain external momenta by using momentum conservation for a given topology, for instance This modifies the options for the function FormCalc CalcFeynAmp so that for LoopNumber = 1 and topology penguinT the momenta of the second external particle is replaced by the momenta conservation expression.
The setting mass allows to neglect masses, and improve the readability of the code, if desired: The first example appends an additional replacement rule that allows the use of the explicit expression for the mass of the particle in the generic propagator.The second one prevents some simplification which would otherwise lead to incorrect amplitude expressions, see the explicit implementation for more details.
Finally, regularization overrides the option FormCalc Dimension for given topologies: This option might be useful to obtain the desired Wilson coefficients faster or more optimally: e.g.sometimes the option chains might be skipped, as the default FormCalc setting has already produced the required expressions.

Example 3: lepton self-energy with NPointFunctions (Observable -> O i [])
Let us exemplify how one enables the usage of the advanced settings for NPointFunctions, by modifying Example 3 which calculates fermion self energies.We continue from Section 3.3.3and Section 3.4.3,where the self-energy calculation was defined by using only the simple mode of NPointFunctions.In order to use the advanced settings, at first we need to apply the following modifications in the file FlexibleSUSY.m:Here the first lines are as explained in Section 3.5.1 and enable the advanced settings, including the variable contr which allows users later to select different Feynman diagrams for the evaluation of the self energies via keywords corresponding to different Contribution i .In the last line of the listing, the name of the function is changed to reflect different possible values of the variable contr.The goal in this example is to allow users to compute "self energies" by including either only one-particle irreducible diagrams, or only diagrams with a tadpole part, or both.Hence we intend to use the advanced settings to define three keywords and associate them with the appropriate topologies via the topologies setting described in Table 6.
In the following we describe a typical workflow how one might interactively obtain the relevant information to create the advanced settings file meta/Observables/O i /NPointFunctions.m which achieves this.We begin with the definition of topologies we want to enable/disable for the calculation.
As the self-energy process is 1 → 1, we start by looking for already defined topologies inside meta/NPointFunctions/Topologies.m in the form AllTopologies[{1, 1}].Currently, this T3: definition is missing, so its specification becomes our first task.We can open a Mathematica notebook and evaluate the following code to figure out, which one-loop 1 → 1 topologies exist in general and which ones are of interest to us: Listing 27: Execute in a Mathematica notebook.
The printed output is similar to the content of the first column in Figure 3.In this example, we are not interested in the first topology, while others represent what would like to include.We need to convert the FeynArts representation of the topology to the NPointFunctions one, this can be done as shown in Figure 3 or by execution of the code: 13Listing 28: Execute in a Mathematica notebook.

AdjaceTopology[topologies[[2]]] AdjaceTopology[topologies[[3]]]
Let us name the obtained topologies as in Figure 3 and add their definition to the observableindependent file meta/NPointFunctions/Topologies.m as: where the last line creates the synonym for both topologies combined.Now we can finally define the content of the advanced settings file for this observable: This achieves the goal.Now the observable ExampleLeptonSE has an option contr, like all the examples in Section 2. This option may be set to the values Tadpoles, Sunsets, or Fermi.If some model-specific configuration file models/M a /FlexibleSUSY.m contains the option Tadpoles, then only the tadpole topology will be included, and similarly for Sunsets.Once we specify Fermi (or {Tadpoles, Sunsets}), then both topologies will be used to compute the (generalized) selfenergy.In this setup, Sunsets has the same effect as the Irreducible setting from the simplified example.

Content of the optional file O i /FSMathLink.m
The spectrum generators created by FlexibleSUSY can be called from within a Wolfram Language notebook or kernel.The necessary definitions to output the numerical values of the generated observables at the Wolfram Language level are done in the general file meta/FSMathLink.m.For a specific observable O i one can create an optional file meta/Observables/O i /FSMathLink.m to specify the desired interface via the function PutObservable.We refer the reader to existing example files shipped with FlexibleSUSY for more details.

General case
To activate the C++ code generation for a desired observable with FlexibleSUSY, one carries out the same steps as for predefined observables in Section 2 (though now it is clear where all definitions come from).For a user-selected model M a one modifies models/M a /FlexibleSUSY.m by changing ExtraSLHAOutputBlocks.The pattern of observable O i is defined inside the file called meta/Observables/O i /Observables.mwhile the Les Houches blocks where this observable can be placed -inside meta/Observables/O i /WriteOut.m.Finally, by execution of make one obtains the desired C++ spectrum generator.

Example 1: output a given number
In summary, to add the new observable corresponding to Example 1 to FlexibleSUSY, one needs to create several files: After all this is done, the observable can be used in the same way as the predefined observables discussed in Section 2, and we need to proceed as described in that section.
We first need to choose a desired physical model to perform the calculations (to continue this example we choose the SM), create it via ./createmodel--name=SM, then configure FlexibleSUSY to make the C++ spectrum generator for this model via ./configure--with-models=SM.
Then, two model-specific settings files have to be modified to enable the computation of the observable ExampleConstantObservable. To be specific, we modify the file with the meta-level model settings models/SM/FlexibleSUSY.m to contain In this way we specify that the observable is called twice, with two different arguments (the arguments simply correspond to the numeric constants which should be printed in the output), and that the output will be part of the Block FlexibleSUSYLowEnergy with numbers 1 and 2, respectively.Finally, the calculation of all observables is enabled in the runtime model-specific SLHA input file: The execution of make command runs the meta phase and compiles the final C++ spectrum generator.FlexibleSUSY provides shell scripts which merge these steps, as follows: Listing 33: Execute in terminal from the FlexibleSUSY directory.
. Here we see the appearance of the block numbers and the values of the numeric constants in agreement with the definitions given above.

Example 2: show fermion masses
In order to add the observable of Example 2 to FlexibleSUSY and obtain a SM spectrum generator including the output of this observable one needs to go through analogous steps.Again, FlexibleSUSY is shipped with a script which merges all steps: Listing 34: Execute in terminal from the FlexibleSUSY directory.
./createmodel --name=SM ./configure--with-models=SM ./examples/new-observable/make-observableSM-example-2 make ./models/SM/run_ SM.x --slha-input-file=models/SM/LesHouches.in.SM After executing these steps successfully, the SLHA output will contain the results corresponding to four instances of the observable ExampleLeptonMass (where the fermion index 1 here means 2nd generation due to the C++ numbering convention): Block FlexibleSUSYLowEnergy Q= 1.73340000E+02 1 1.04187667E-01 # Fe [1] (lepton [1] if in Block ExampleLeptonMass) mass 2 5.71258815E-02 # Fd [1] (lepton [1] if in Block ExampleLeptonMass) mass Block ExampleLeptonMass Q= 1.73340000E+02 1 1.05658357E-01 # Fe [1] (lepton [1] if in Block ExampleLeptonMass) mass 2 1.05658357E-01 # Fd [1] (lepton [1] if in Block ExampleLeptonMass) mass Here the file models/SM/FlexibleSUSY.m defines the four instances of the observable: they are distinguished by their argument (either Fe [2] or Fd [2]) and by the Les Houches block (either FlexibleSUSYLowEnergy or ExampleLeptonMass).We recall that in the definition of the observable the running MS/DR mass of the fermion in the BSM model described by the argument and the SM lepton pole mass of the specified generation (in this case generation 2) are returned, see Section 3.4.2.We also recall that the output depends on the Les Houches block, see Section 3.2.2.This explains the output given above, where the FlexibleSUSYLowEnergy block shows the MS/DR masses of the muon and the strange quark in the BSM model, while the ExampleLeptonMass block shows twice the muon pole mass from the qedqcd object.

Example 3: lepton self-energy with NPointFunctions
Technically, there are two versions of this example which differ by O i /FlexibleSUSY.m: using the simple mode of NPointFunctions in Section 3.3.3or using the advanced mode in Section 3.5.7.To add the advanced version of the observable of Example 3 to FlexibleSUSY and obtain a SM spectrum generator, one needs to go through steps analogous to previous examples.Again, FlexibleSUSY provides a script which combines all steps: Listing 37: Execute in terminal from the FlexibleSUSY directory.

Conclusions
In this paper, we describe two novel essential features of FlexibleSUSY: a streamlined approach for incorporating new observables to FlexibleSUSY, and a simplified way to generate C++ code for Feynman diagrams essential for the computation of these observables and other physical quantities.To enhance the accessibility of this article we have divided the content into two distinct parts.
The first part is concise, requires no detailed knowledge of the code and is directed to readers interested in the computation of predefined FlexibleSUSY observables.The currently predefined observables include CLFV observables such as ℓ i → ℓ j γ, µ-e conversion and ℓ i → ℓ j ℓ k ℓ c k , whose physics definitions are summarized in the appendix.Any user of FlexibleSUSY can now enable the computation of these observables for any desired model by just adding the appropriate flag in the model files.
The second part deals with the procedure of implementing new observables.It is more extensive and technical and is of interest to those seeking insights into the internal structure of the NPointFunctions code used for the creation of new observables.In essence, the implementation of a new observable requires five files (three Wolfram Language and two C++ template files) which essentially define in a model-independent way how the observable is computed, how it is called and how its output is organized.The second part also illustrates the utilization of the NPointFunctions extension to streamline the generation of C++ code for Feynman diagrams.This is achieved with the help of FlexibleSUSY-specialized wrappers designed for FeynArts, FormCalc, and ColorMath packages.
Three fully worked out examples are provided.They correspond to a minimal "observable" which simply outputs a single number, an observable which outputs a set of particle masses, and an observable which outputs the values of one-loop self energy diagrams.The code snippets can be used as efficient starting points for future implementations of further observables.In the appendix additional features and details are discussed which may be useful to accommodate special properties of new observables.They correspond to actual use-cases motivated by implementing e.g. the b → sµµ or h → gg decays.
The FlexibleSUSY extensions presented here have been thoroughly cross-checked and used for concrete phenomenological applications in non-SUSY and in SUSY models.Ref. [23] uses and validates CLFV observables in a leptoquark model, where some observables arise at the one-loop level and some observables arise at tree-level.Refs.[24,25] provide validations of loop-induced CLFV observables in a model of neutrino masses where several neutrino masses are themselves loop-induced.Finally, extensive validations have been carried out in the context of a non-trivial SUSY realization by comparing with results of Ref. [26] on CLFV phenomenology in the MRSSM.These applications demonstrate the reliability and versatility of the code for observables defined via Feynman diagrammatic calculations across a broad spectrum of models for wide variety of observables.In the future, it is planned to add further observables to the default distribution of FlexibleSUSY.In addition, users may add observables individually.Finally, we introduce the possibility to request the implementation of desired observables via github issues, see the developer's repository.
The observable defined in this example calculates amplitudes required for the h → gg process (but not the branching ratio itself) with the help of NPointFunctions.After NPointFunctions calculations are done, one is left with amplitudes containing many abbreviations.From physics we know that the computed amplitudes must contain several structures of covariants.For h → gg, the possible covariants are ǫ µ 2 ǫ 3µ , ǫ µ 2 p 3µ , ǫ µ 3 p 2µ , ε αβµν ǫ 2α ǫ 3β p 2µ p 3ν , (B.1) where ǫ i is the polarization vector of the corresponding gluon and ε is the Levi-Civita tensor.
It is the coefficients of these structures which are relevant for the calculation of the branching ratio.Hence we need to extract the prefactors of these structures.This can be performed as shown in Listing 40.There, we apply all remaining sub-expressions in line 1, then abbreviate all basis structures in lines 2-8, and extract the coefficients that come as prefactors of the second argument of InterfaceToMatching in line 9 (similarly to line 28 in Listing 18).Finally, the C++ code definitions are generated by the function NPFDefinitions in line 14 (similarly to line 32 in Listing 18; note, that NPFDefinitions may accept strings that are different to the last argument of InterfaceToMatching): npf = WilsonCoeffs InterfaceToMatching[npf, {"eps", "e2e3", "e2m3" "e3m2"}]; ...

AppendTo[npfDefinitions,
NPointFunctions NPFDefinitions[npf, "cpp _ name", SARAH Delta, {"eps", "e2e3", "e2m3 _ e3m2"}] ]; In principle, one can apply the routines from the code snippet above for the processes with external fermions as well.This might be relevant, in particular, if one prefers to ignore the setting chains from Section 3.5.5 and deal with Dirac chains in some other ways.

Table 1 :
All observables currently supported by FlexibleSUSY.
lep symbol Leptons ℓ in ℓ i → ℓ j ℓ k ℓ c k in SARAH model (Fe in SM or MRSSM) iI, iJ, iK integer Generations i, j, k in ℓ i → ℓ j ℓ k ℓ c k , starting from 1

2
Listing 2: In models/Ma /LesHouches.in.Ma .All added observables are enabled.FlexibleSUSY One can also output Wilson coefficients used in derivation of LToLConversion or BrLTo3L.To do that, one places the corresponding observable into the FWCOEF (IMFWCOEF) block to output their real (imaginary) part, for example:
{ "num _ value" -> "Re(observables."<> Observables GetObservableName[obs] <> ")", "leptons" -> Switch[gen, 0, "1111", 1, "1313", 2, "1515"] } ]; 1as indicated in the listing above).In general, we are interested in defining one observable O i unified by similar calculations, e.g.we define one observable O i for the set of processes ℓ i → ℓ j ℓ k ℓ c k instead of multiple observables µ → 3e, τ → 3µ, etc.Then, we enable specific realizations of O i in ExtraSLHAOutputBlocks, see also explicit examples in Listing 31 and 35.So, we start with selecting unique realizations of chosen O i from all allObs and storing them into observables.Then, inside If statement we make changes to the C++ code generation based on the possible realization-specific features (e.g. the process τ → 3µ might require additional contributions, compared to µ → 3e).Changes in C++ prototypes stored in the variable prototypes are handled automatically due to the function WriteObservable, see Listing 5, while the C++ definitions in the variable definitions usually require manual coding, see Listings16-18.

Table 4 :
Mandatory options for the function NPointFunction, see the function CheckOptionValues in meta/NPointFunctions.m file for allowed values, and meta/Observables/Oi /FlexibleSUSY.m files for examples.
3.3.For this reason we use the C++ template files given in

Table 5 :
Purpose of all available advanced settings and their usage.any changes in this example.Note that these general template files are generally appropriate for observables that use the NPointFunctions extension, hence we can keep the present subsection very short.

Table 6 :
Table 6 shows the required command in the file O i /NPointFunctions.m.It connects each keyword Contribution i with certain Syntax and options values for topologies.

Table 6 WhenToApply
treeS, and NPointFunctions internal representation of this specific topology.It is obtained in a few steps, as follows Present or Absent Apply Command b for TopologyName j if Contribution i is present (absent) in KeepProcesses, see Table 4

Table 7 :
Syntax and options for diagrams and amplitudes.