SpecBit, DecayBit and PrecisionBit: GAMBIT modules for computing mass spectra, particle decay rates and precision observables

We present the GAMBIT modules SpecBit, DecayBit and PrecisionBit. Together they provide a new framework for linking publicly available spectrum generators, decay codes and other precision observable calculations in a physically and statistically consistent manner. This allows users to automatically run various combinations of existing codes as if they are a single package. The modular design allows software packages fulfilling the same role to be exchanged freely at runtime, with the results presented in a common format that can be easily passed to downstream dark matter, collider and flavour codes. These modules constitute an essential part of the broader GAMBIT framework, a major new software package for performing global fits. In this paper we present the observable calculations, data, and likelihood functions implemented in the three modules, as well as the conventions and assumptions used in interfacing them with external codes. We also present 3-BIT-HIT, a command-line utility for computing mass spectra, couplings, decays and precision observables in the MSSM, which shows how the three modules can be easily used independently of GAMBIT.


Introduction
Run II of the Large Hadron Collider (LHC) is engaged in a wide-ranging search for evidence of physics Beyond the Standard Model (BSM). Such models typically have large parameter spaces, so understanding their phenomenology requires detailed calculations using computer programs. Phenomenological studies therefore often involve a large chain of public codes that must be linked together.
This set of codes includes spectrum generators (to determine the masses of new particles), decay calculators (to obtain decay widths), and packages capable of predicting low-energy precision observables, such as the anomalous magnetic moment of the muon. These codes need to be linked together in such a way that information from the spectrum generator can be passed to the other calculators, and their outputs can in turn be used in other programs.
For the Minimal Supersymmetric Standard Model (MSSM) and the Next-to-minimal Supersymmetric Standard Model (NMSSM), there are the SLHA [1] and SLHA2 [2] conventions, which simplify matters somewhat. However, even in these cases it can be far more convenient to have the programs automatically linked, as can be testified by the popularity of packages that incorporate individual codes, such as SUSY-HIT [3] for the MSSM and NMSSMTools [4][5][6] for the NMSSM. On the other hand, keeping individual codes distinct in a modular framework allows for easier isolation of the origin of differences in results, and for the ability to mix and match different codes according to preference or the specific advantages of one tool or another. Indeed, the need to compare several software packages for cross-checks, catching bugs and revealing uncertainties in each calculation has been demonstrated on numerous occasions (e.g. Refs. [7][8][9]).
Here we present a framework that provides the best of both worlds, designed to work with all models including both non-supersymmetric models and exotic SUSY models beyond the MSSM and NMSSM. The framework consists of three packages: SpecBit, for handing renormalisation group running and calculation of mass spectra, DecayBit, for computing branching ratios and widths, and PrecisionBit, for calculating other precision observables. All three packages are interfaced in a common way within GAMBIT 1.0.0 [10], allowing specific functions and external software to be exchanged at runtime, whilst being run from exactly the same input parameters under the same physical conventions. The user interface is designed to be as simple and general as possible, so that information can be extracted and inserted in an intuitive way.
GAMBIT (the Global and Modular Beyond-SM Inference Tool) is a multi-purpose physics tool for performing parameter scans and global fits in BSM models using either frequentist or Bayesian statistics. Each of the three packages we describe here is a module within the broader GAMBIT framework. Here we highlight terms that may be considered GAMBIT jargon, and provide a simple summary of their meanings in the glossary (Appendix H). Each module collects a series of module functions based around a common theme.
Each module function typically calculates a single observable or likelihood component for use in a fit, or performs some calculation that is needed by another module function in order to eventually arrive at the value of an observable or likelihood. Module functions that require the results of other module functions can declare dependencies on the physical or mathematical quantities that they require. GAMBIT defines a series of theoretical models available for analysis, and various rules that relate the different module functions to each other, to models, and to functions that are available for modules to call from external physics codes, known as backends. At runtime, the GAMBIT dependency resolver identifies the module and backend functions required to compute the observables and likelihoods requested by a user, and arranges them into a dependency tree. It then uses methods from graph theory to 'solve' the tree, and determine the order in which the functions must be called so that all dependencies and model consistency requirements are guaranteed to be satisfied. GAMBIT's sampling module ScannerBit [11] runs the user's choice of sampling algorithm on the solved dependency tree, and saves the resulting parameter samples, derived observables and likelihoods in an output database for subsequent statistical analysis and plotting.
SpecBit, DecayBit and PrecisionBit are used to compute and return spectra, decay widths and precision observables, run associated backend codes, interpret the results thus obtained, and provide them to other modules as required for subsequent calculations. Other GAM-BIT modules calculate dark matter observables (DarkBit [12]), high-energy collider signatures (ColliderBit [13]) and quantities from flavour physics (FlavBit [14]). An extended description of the structure, features and abilities of GAMBIT can be found in Ref. [10].
SpecBit, DecayBit and PrecisionBit are designed so that a user may easily add new models and interfaces to new backends. Models that have already been implemented in GAMBIT are automatically supported. The first release features various incarnations of the MSSM and the scalar singlet dark matter model, as these are the models on which the first GAMBIT scans [15][16][17] have been performed.
In Sec. 2 we describe the module SpecBit. This includes how SpecBit manages the running of spectrum generators, how to access and use the information it extracts from them, and how it can be used to analyse vacuum stability. Next we provide a detailed description of the module DecayBit (Sec. 3), including the decay data it contains, decay calculations it performs internally, and its use of backends. In Sec. 4 we detail Precision-Bit, including its likelihood functions and interfaces to external precision codes. Sec. 5 gives some examples of each of the modules in action within GAMBIT, and presents 3-BIT-HIT, a simple utility that performs a similar function to SUSY-HIT, and serves as a basic example of the standalone use of SpecBit, DecayBit and PrecisionBit outside of the GAMBIT framework. We give a brief summary in Sec. 6. In Appendix A we review the physics of the models discussed in this paper, and the conventions that we adopt for them. We provide details of the interface with SPheno in Appendix B, explicit documentation of some of the classes involved in SpecBit in Appendices C-E, instructions for adding new spectrum generators to SpecBit in Appendix F, and an explicit example of how to build spectrum classes for new models in Appendix G. In Appendix H, we give the glossary of GAMBIT terms highlighted at various points in this paper.
SpecBit, DecayBit and PrecisionBit are released under the standard 3-clause BSD license 2 , and can be downloaded from gambit.hepforge.org.

SpecBit
All information in GAMBIT about the spectrum of particle masses and their couplings comes from module Model parameters are obtained from the GAMBIT Core and fed to SpecBit module functions. These module functions run spectrum generator backend codes and embed the results within Spectrum wrapper objects, which carry the spectrum information out to other GAMBIT functions that require it. In some cases the backend code remains connected to the wrapper object, allowing RGEs to be called to translate couplings and mass parameters to other scales. In this diagram solid outlines indicate "active" elements for a hypothetical scan, while dashed outlines indicate inactive elements. The "active" elements are those that are activated by the central GAMBIT dependency resolution system, or manually run in a standalone code. The "inactive" elements represent alternate calculation pathways not required in a scan, but can be switched on instead if the user chooses.
functions in SpecBit. This includes the pole masses and mixings of all physical states in the model, DR [47][48][49] or M S [50] parameters, Higgs couplings and basic SM inputs. In BSM models where some of this required information is a calculable prediction of the model, SpecBit will obtain it by calling a spectrum generator, taking inputs as model parameters. In cases where the information is not a prediction, but is already specified directly by the free parameters of the model (e.g. if a pole mass is defined as a model parameter), SpecBit will simply take this information from the model parameters and store it in a similar way to the information extracted from the spectrum generator. Spectrum data is then stored inside Spectrum wrapper objects, which transport the data to other GAMBIT module functions that request it. A general overview of this process is shown in Figure 1.
Because some calculations require the values of running parameters at a particular scale, SpecBit also allows the DR or M S parameters to be run to a scale chosen locally in any module function, using a relevant spectrum generator backend.

Supported models and spectrum generators
The first release of SpecBit is extendable to any model, but has built-in support for the scalar singlet dark matter model, and the MSSM. It also provides a low-energy spectrum object containing SM information.

Standard model spectrum
GAMBIT is designed to fit models of BSM physics. There is thus no full spectrum generator implemented for the SM. However, BSM spectrum generators typically rely on SM inputs as low energy boundary conditions. To store this information, SpecBit contains a QedQcd spectrum object, which was originally part of SOFTSUSY and is also used in FlexibleSUSY. This is the source of some of the SM data that can be extracted from SpecBit 3 and details of the relevant capabilities are given in Sec. 2.2.2 and Table 1. More details on this QedQcd object can be found in the SOFTSUSY manual [30].

Spectrum generators for the MSSM
MSSM mass spectra are typically obtained by finding solutions to the renormalisation group equations (RGEs) that simultaneously satisfy boundary conditions (BCs) at high and low scales, before using self energies to calculate the pole masses. The low-scale BC matches the spectrum to observed SM data, and the high-scale BC places constraints on the soft SUSY-breaking masses.
The first version of SpecBit comes with interfaces to two different spectrum generators: FlexibleSUSY and SPheno. It also contains an interface to FeynHiggs, which can be used to obtain Higgs and sparticle pole masses; details of our interface to FeynHiggs are given in Sec. 3.1.3. FlexibleSUSY creates a spectrum generator for a given model, defined in a SARAH input file, and uses a FlexibleSUSY model file to specify the BCs for that model. SPheno runs in different modes according to the input BCs.
The model is defined by both a set of input parameters that specify DR parameters of the MSSM, and an additional scale at which the soft SUSY-breaking DR parameters are defined. Note that when there are constraints relating the parameters to each other, varying this scale is not equivalent to a reparameterisation, as there will be different mass splittings that cannot be reproduced by any point in parameter space when the constraints are applied at a different scale.
In the spectrum generator MSSMatMGUT, the soft parameters are defined at the scale where the gauge couplings unify, which is determined iteratively. In MSSMatQ, the soft parameters are instead defined at a user-specified scale Q. The MSSMatMGUT and MSSMatQ spectrum generators support the MSSM63atMGUT and MSSM63atQ models respectively, which are the most general options in the GAMBIT model hierarchy described in Sec. 5 of Ref. [10].
More constraining BCs create lower-dimensional subspaces of these more general parameter spaces, so a vast number of possible subspaces exist. SpecBit uses the GAMBIT model hierarchy (see Sec. 5 of Ref. [10]) to relate these subspaces to the MSSM63atMGUT or MSSM63atQ models, and can use the MSSMatMGUT or MSSMatQ spectrum generators to determine their mass spectra. It is also possible to directly implement the boundary condition in a FlexibleSUSY spectrum generator. As an example of this, we have also implemented a specific CMSSM spectrum generator. This was mostly introduced as a basic check of the the model hierarchy, but remains in GAMBIT for convenience.
Technical details of how the FlexibleSUSY spectrum generators are implemented are given in documentation shipped with the code (see doc/Adding_FlexibleSUSY_Models.txt), though an illustrative example is given in Appendix G.1, and details required to use them in scans are given in sections 2.2.4 and 2.2.5. In the additional documentation, we also demonstrate how to add an MSSM variant where soft parameters are fixed at the SUSY scale (lowMSSM in the FlexibleSUSY naming scheme), where the SUSY scale is defined as the geometric mean of the DR stop masses and is determined iteratively. We will include this variant natively in the next version of SpecBit.
In the case of SPheno, the running mode is determined by the initialisation of its internal variables. In contrast to FlexibleSUSY, the available modes are triggered according to the model being scanned (CMSSM, or one of the MSSMatMGUT or MSSMatQ models), not by the usage of different functions. Other running modes available in the out-of-the-box version of SPheno, such as the NMSSM, GMSB, AMSB, etc., are not covered in the backend version of the software used in GAMBIT but, as mentioned before, there are immediate plans to include these in future releases. The specific details of how to use SPheno as a spectrum generator in GAMBIT are provided in Secs. 2.2.4 and 2.2.6. A comprehensive description of the backend system can be found in the main GAMBIT paper [10], and in Appendix B for the specific case of SPheno.
For each of these spectrum generators the EWSB conditions are used to fix |µ| and b, so that the Higgs vacuum expectation values (VEVs) are fixed by the measured Z boson mass and the input ratio of the two VEVs (tan β = vu v d ).

Spectrum generators for the scalar singlet dark matter model
The scalar singlet mass spectrum is available in two forms. One is a simple container object that is setup without any spectrum generation, where the parameters may be interpreted as being either M S parameters, in which case the pole masses will only be known at tree level, or they could be interpreted in the OS scheme. The other option is a spectrum object that interfaces with FlexibleSUSY to calculate the spectrum.
The simple container spectrum is the most efficient option for a range of calculations in GAMBIT that require only the masses and couplings at a fixed scale. In this case, the only new model parameters required for other module functions, m S and λ hs , are free to vary independently of the other SM parameters, as any effect on the SM contents only occur at the one-loop level. In the container object, the pole masses and couplings are stored directly from the GAMBIT model parameters. Other quantities not available in the model parameters, such as the gauge and Yukawa couplings, are calculated from SM parameters at tree level.
If radiative corrections and/or renormalisation group running of model parameters is necessary, then a fullycalculable spectrum is required. For this purpose, we offer a spectrum object that uses FlexibleSUSY and the input parameters to calculate pole masses and couplings. This spectrum object also has the capacity to evolve parameters between different scales using RGEs. This spectrum object interprets the input model parameters, in particular the masses of the scalar singlet (m S ) and the Higgs boson (m h ), as M S quantities, while all other SM masses are taken according to SLHA2 conventions [2]. These running parameters are defined at the scale m Z . The EWSB conditions are then imposed to calculate the value of the quartic Higgs coupling λ (Eq. A.2), and pole masses are calculated at a user-specified scale Q. The scale Q may be varied as a nuisance parameter in the GAMBIT model SM_Higgs_running (see Ref. [10] for more details), or left as the top mass by default. Optimally, this scale should be set to minimise logarithmic contributions to the RGE β functions with the largest coefficient.

General settings
In this section we describe how to use SpecBit to run available spectrum generators via a GAMBIT YAML 4 file. At the most rudimentary usage level, spectrum information can simply be written to disk for one model point in an SLHA2-like format, using SLHAea 5 (a C++ class allowing internal representation of SLHA files). At a more advanced level, the information can be written to disk for every point in a GAMBIT scan, via the GAMBIT printer system (see Ref. [10]), and analyzed using external software. More advanced usage, such as accessing spectrum information at the C++ level (in a GAMBIT module function, for example), is described in section 2.3. For details on YAML usage of the capabilities and associated options described here, please see the GAMBIT manual [10] or the SpecBit example files in the yaml_files directory of the GAMBIT source tree. A README file can be found in this directory, which explains the example files further.

Standard model
The capabilities available in SpecBit related to the SM are given in Table 1. The capability SM_spectrum is provided by the function get_SMINPUTS.
This creates an object containing the low-energy SM parameters obtained directly from the model parameters, which must contain the parameters in the StandardModel_SLHA2 GAMBIT model [10]. This can then be used to build BSM spectrum objects that require this low-energy data.
By default, the get_SMINPUTS function populates its W mass with the observed value of 80.385 GeV [51]. Users who prefer the W mass returned to respect the tree-level relationship with m Z and sin 2 θ W can choose for it to be calculated at tree level instead, by setting the option enforce_tree_level_MW = true for this function.
The module function get_QedQcd_spectrum has a dependency on SM_spectrum, and provides the capability qedqcd_subspectrum. This function creates an effective QedQcd object. The returned QedQcd object contains the low energy data, where all running quantities from the SMinputs dependency are now given at the scale m Z . This is then wrapped along with the original SMinputs into a SubSpectrum object. 4 www.yaml.org 5 fthomas.github.io/slhaea Finally, to provide the capability SM_spectrum, the module function get_SM_spectrum wraps the qedqcd_subspectrum object into a full spectrum object along with a simple container for the Higgs pole mass and vacuum expectation value, providing a complete SM spectrum.

Scalar singlet dark matter
The capabilities available in SpecBit that either produce or depend on the scalar singlet spectrum are given in Table 3.
SpecBit has a capability SingletDM_spectrum to provide essential details about the scalar singlet spectrum. This capability is provided by the module functions get_SingletDM_spectrum_simple and get_SingletDM_spectrum_FS, both of which return a Spectrum object. The former returns the simple container spectrum, as described in Sec. F.2, and the latter returns a more sophisticated spectrum, computed with FlexibleSUSY, which can run parameters as described in Sec. 2.3.3. All options described in Sec. 2.2.5 can be given to this function.
The validity of the spectrum is checked by Flex-ibleSUSY during the initial calculation, and may be optionally checked up to arbitrarily high energies afterwards. If an error is encountered during initial spectrum generation, a FlexibleSUSY error is passed to GAMBIT. The handling of such an error is controlled by the option invalid_point_fatal; if this option is set to true, invalid specctra will trigger a SpecBit error. It is a common consideration to check for perturbativity of the dimensionless couplings up to a specific high energy scale [52][53][54]. We provide the ability to do this with the option check_perturb. The maximum scale to run the couplings up to is set with the option check_high_scale. If the couplings are found to be non-perturbative, then an invalid point exception is raised as if the initial spectrum calculation encountered an error. The options, input types and default values for the get_SingletDM_spectrum_FS function are given in the list below.
-check_perturb: takes a bool to demand that the spectrum be run to check_high_scale. Default false. The module function get_SingletDM_spectrum_as_map also provides capability SingletDM_spectrum, but as a C++map (aliased as map_str_dbl). For this function to run, the SingletDM_spectrum needs to be provided as a Spectrum object. In other words, this function just translates between the Spectrum type and a C++ map. The main use of this function is to print the contents of the spectrum to an output stream for each data point during a scan (as the map_str_dbl type is "printable" but the Spectrum type isn't), as shown in Fig. 2.
There are also some additional capabilities related to vacuum stability in the scalar singlet model: vacuum_stability, check_perturb_min_lambda, VS_likelihood and expected_lifetime. These are discussed in detail in the advanced usage example in Sec. 2.5.2.

MSSM
The capabilities available in SpecBit relevant for the MSSM model are given in Table 2.
SpecBit has a number of module functions that provide the capability unimproved_MSSM_spectrum. We will start by examining three of them:  also provides this unimproved_MSSM_spectrum capability in the form a Spectrum object. However, instead of using a spectrum generator to do this it uses an SLHA file. The Spectrum object in this case does not provide the ability for RGE evolution, but is rather a simple container spectrum. This function takes the options filenames: path to SLHA file to be read in. Can be a list of many file names, in which case each file will be read in sequentially. cycles: number of loops over the filenames list to allow (default -1 indicates no limit). When the limit is hit an error will be raised, stopping the scan.
This function is mainly for debugging purposes, because one usually generates spectra on the fly during a scan. Two further SpecBit module functions, get_unimproved_MSSM_spectrum_as_SLHAea and get_unimproved_MSSM_spectrum_as_map are available for translating a Spectrum object to either an SLHAea object (SLHAstruct) or to a C++ map (as in the analogous scalar singlet DM function), respectively. These still represent the spectrum, just with a different type, so the module functions have both a dependency and a capability of unimproved_MSSM_spectrum, but with different types. See Fig. 2 for an example use case of this 'translation' behaviour. We have found each format useful for different tasks. For example, some module functions in other parts of GAMBIT request spectrum information in SLHAstruct format, rather than as full Spectrum objects, because these functions work closely with backend codes that were originally designed to work via SLHA files. So keeping the format directly parallel to SLHA was convenient for them. On the other hand, the std::map format is useful because at present neither Spectrum objects nor SLHAstruct objects can be directly written to disk by the GAMBIT printer system [10], so converting to std::map is necessary as an intermediate step for this purpose. The module functions get_MSSM_spectrum_as_SLHAea and get_MSSM_spectrum_as_map also translate a Spectrum object into an SLHAea structure or a std::map. However in these cases the capability and dependency is MSSM_spectrum, a capability that can also be provided by the module functions from PrecisionBit, as shown in module function Table 10 and discussed in Sec. 4.2.3.
In the MSSM it is often useful to be able to identify which Higgs boson is most similar to the SM Higgs boson. The module function most_SMlike_Higgs_MSSM provides this information by returning a PDG code representing the CP -even Higgs with couplings that are closest to those of the SM Higgs boson. We deem the lightest CP -even state to be the most SM-like if sin(β − α) > cos(β − α), where α is the mixing angle defined in Eq. A.16; otherwise, this function returns the PDG code of the heavier CP -even Higgs. This function has capability SMlike_Higgs_PDG_code, and a dependency on the MSSM_spectrum.
It is also possible to extract the low-energy SM spectrum from the MSSM spectrum. The capability SM_subspectrum is provided by the module function get_SM_SubSpectrum_from_MSSM_Spectrum, which has a dependency on the MSSM_spectrum.
Finally, there are three module functions that require backend functions from FeynHiggs. These module functions,

FlexibleSUSY options
Any spectrum generator interfaced to GAMBIT will inevitably have its own set of options to control precision and methods of calculation. Here we briefly summarise the FlexibleSUSY options, which can be set via YAML options. For a more detailed explanation of these options, see the FlexibleSUSY manual [18].
x /(4π) and y x is the corresponding Yukawa coupling, and α s is the strong coupling constant.
If the FlexibleSUSY spectrum generator encounters an error during calculation, this may be passed on as a SpecBit error, resulting in the termination of a scan, or as an invalid point exception, which will result in the scan point being given an invalid likelihood (an extremely small value defined by model_invalid_for_lnlike_below in the YAML file). See Ref. [10] for more details on exception handling. The management of these exceptions within GAMBIT is controlled with the option invalid_point_fatal discussed in Sec. 2.2.3.

SPheno options
As with FlexibleSUSY, there are a number of options available for SPheno to control certain aspects of the calculation. In GAMBIT only a limited set of all the available options for SPheno is allowed, because many features or models permitted by SPheno are not yet covered in GAMBIT. A detailed explanation of all the options available for SPheno can be found in the manual [20], under the SLHA block SPhenoInput.
The following options can be provided via the YAML options section in the input YAML file, and they will be internally assigned to the corresponding variables in SPheno. Errors prompted by SPheno are collected by the internal integer variable kont, which takes different values according to the source of the error. The specific error messages corresponding to values of kont can be seen in Appendix C of the manual [20]. For any of these values, the backend convenience function run_SPheno raises an invalid_point exception, as described in [10].

Mass cut options
For any module function that constructs a Spectrum object, it is possible to specify options that enforce user-defined relationships between particle pole masses: mass_cut and mass_cut_ratio. Any spectrum that does not pass the specified cuts is declared invalid (which may invalidate the entire parameter point, depending on whether the likelihood ultimately depends on the spectrum or not [10]). While these cuts are not, in general, physically required, we have found this feature useful for certain specialised scans, for example explorations of co-annihilation regions in the MSSM. These options take one or two particle names recognised by the spectrum, along with two numbers used to define the range to be cut. They are best explained by way of an example entry in the Rules section of a GAMBIT YAML file: The above code invalidates any points that do not have: the mass of the lightest Higgs between 120 and 130 GeV, -the mass of the lightest slepton between 100 and 1000 GeV, and the mass of the lightest slepton within 1% of the mass of the lightest neutralino.
Note in particular the use of absolute value signs for the neutralino mass, which may be negative. This notation can be used wherever a particle name is given, in order to define the cut rule in terms of the absolute value of the particle mass. Also note that two such cut rules are given for mass_cut, illustrating the fact that these options can be used to apply arbitrarily many cuts simultaneously.

Interface details for GAMBIT module writers (C++ API for Spectrum and related classes)
In Sec. 2.2 we provided details for a base-level user of GAMBIT to simply operate the code as written. However, one of the goals of GAMBIT is to provide a framework into which researchers can add their own calculations while maintaining easy use of all pre-existing GAMBIT functions. This section is intended as a guide for users of this kind. This information will also be needed by users running SpecBit externally to GAMBIT. Here we describe details of the interfaces to objects provided by the SpecBit module, and how they should be used to facilitate physics calculations of other modules. In Sec. 2.3.1 we demonstrate the most common methods for accessing spectrum information. In Secs. 2.3.2 and 2.3.3, we outline the class structure used to store spectrum information, and explain the details of its interface. The helper class SMInputs is discussed last, in Sec. 2.3.5.

Basic spectrum access
The spectrum information is stored in a Spectrum object, the full structure of which will be described in Sec. 2.3.2. However here we briefly describe the basic user interface.
To access a Spectrum object in a new module function, the module function should be declared to have a dependency on the appropriate spectrum. As described in section 3.2.1 of the GAMBIT paper [10], the spectrum object can then be obtained by dereferencing the safe pointer returned by the SpecBit module function, e.g. using namespace Pipes::my_module_function; const Spectrum& spec = *Dep::SingletDM_spectrum; One can now access the spectrum information using spec. The spectrum information is primarily accessed using a string accessor plus a tag, as in double mh0 = spec.get(Par::Pole_Mass,"h0_1"); Note that Par is simply a namespace, so the tags can be brought into scope with the using keyword, and then used more succinctly, i.e. using namespace Par; double mh0 = spec.get(Pole_Mass,"h0_1"); The tag in the first argument of the getter specifies the type of information to look for. The second argument is a string that, in combination with the tag, tells the getter exactly what information is requested. To simplify access to quantities that are naturally grouped together in vectors (e.g. CP-even Higgs or neutralino masses in the MSSM), the get function may also be called with an index: double mh0 = spec.get(Pole_Mass,"h0",1); As well as a get function, Spectrum objects contain a has function with an almost identical function signature, which can be used to check if a quantity exists in the spectrum before attempting to retrieve it: bool has_mh0_1 = spec.has(Pole_Mass,"h0",1); where has_mh0_1 is true if the quantity exists in the spectrum, and false if it does not.
It is also possible to set the values of spectrum parameters, but not via the base-level Spectrum object interface. For this more advanced usage, see Sec. 2.3.3. Most users will not need to do this.
This form of string getter allows the spectrum information to be accessed in a simple, uniform way across all models and all GAMBIT modules. The user simply needs to know the specific tags and string names used for each piece of information. Particle names, for example, can be found in the GAMBIT particle database (see Sec. 10.1 of Ref. [10]).
For convenience, in Appendix C we give tables listing the tags, strings, and indices needed to access data via the getters, checkers, and setters for the the SM, the MSSM, and the scalar singlet model.
Note that only information tagged with Pole_Mass can be retrieved via the base-level Spectrum interface. For other information (and to run RGEs, when available) one must use the equivalent SubSpectrum interface, which is described in Sec. 2.3.3.

Spectrum class structure
At the centre of the SpecBit module is a virtual interface class named Spectrum, for accessing typical spectrum generator output in a generalised and standardised way. Wherever possible, SLHA2 conventions [1,2] are used as the standard. The objects accessed by this interface can thus be considered as an in-memory representation of the data that one would typically find in the SLHA2 ASCII output of a spectrum generator.
Each Spectrum object contains three main data members, illustrated in Fig. 3: an SMInputs object, which contains SM input parameters (closely mirroring the SMINPUTS block defined by SLHA2); and two SubSpectrum objects, one labelled LE for 'low-energy' and which typically contains low-scale SM information, and the other labelled HE for 'high-energy', which typically contains higher-scale model information such as Higgs-sector or BSM data. The reason for this separation is that, typically, one requires different calculations for accurately evaluating and evolving low-energy SM parameters (like     Standard Model "input" information is contained in a structure call SMInputs (mirroring SLHA2). Low scale (generally < m Z ) spectrum information is wrapped by a "low energy" SubSpectrum object. High scale (generally > m Z ) spectrum information is wrapped in a "high energy" SubSpectrum object. This split mirrors a common requirement of spectrum generators that Standard Model input parameters be first supplied at a common scale, then run by the "low energy" SubSpectrum object to a matching scale, after which the main spectrum generation for the BSM model is performed and wrapped into the "high energy" SubSpectrum object. The outer Spectrum object provides a uniform interface to these underlying structures and makes it easy to retrieve common parameters such as pole masses without the need to interact directly with the subspectra.
quark masses) and high-energy ones like the MSSM. The code to do each calculation will be wrapped in a separate SubSpectrum structure. This wrapper also provides access to running parameters and RGE-running facili-ties, if these are available in a backend code wrapped by a particular SubSpectrum. We describe the SubSpectrum interface in Sec. 2.3.3. Details needed for writing such a class are found in Sec. 2.4. The purpose of the host Spectrum object is to provide a consistent interface to the underlying subspectra and SM input parameters. In Fig. 4 we provide a simplified overview of the contents of the Spectrum class, for reference. We will now discuss the various functions and data members shown in this figure.

Data member getters
The central members of a Spectrum object are accessed via three functions: get_LE()/get_HE() Accesses the hosted SubSpectrum object identified as 'low-energy' or 'high-energy'. get_SMInputs() Accesses the hosted SMInputs object. clone_LE()/clone_HE() Creates a copy of the hosted SubSpectrum object identified as 'low-energy' or 'highenergy'.
The 'clone' functions exist because by design, it is not possible to perform actions that modify contents of Spectrum objects provided to GAMBIT module functions as dependencies. This protection extends to the SubSpectrum objects, so in order to perform actions that modify the spectrum (like RGE running) one must first copy the whole Spectrum or 'clone' the relevant SubSpectrum.

Fig. 4:
Simplified class declaration of primary Spectrum class. Not shown are various overloads, constructors and assignment operators, helper functions, and all private members. Additionally, the index arguments for the has/get functions are implemented via operator overloads rather than optional arguments as is shown. For the full class declaration please see Elements/include/gambit/Elements/spectrum.hpp in the GAMBIT source code.

Copying Spectrum objects
Because GAMBIT prevents the results of module functions from being modified by other module functions, if a user wishes to modify the contents of a SpecBitprovided Spectrum object, the user must first copy it. The copy constructor (and associated constructors) of the Spectrum object are designed to perform a 'deep' copy. To make a copy, one can simply use the copy constructor or assignment operators, i.e.

RGE running
There is little in the underlying object that can be modified from the Spectrum interface, except via the RunBothToScale function. This function runs the renormalisation-scale-dependent parameters of the underlying SubSpectrum objects to the same scale. Note that the underlying RGEs may not be valid beyond certain scales on both the high and low ends, so some caution and knowledge of the underlying objects is required in order to safely use this feature (See 'RGE running' in the following subsection). These objects interface directly to spectrum information or external spectrum generators. Accessor functions (e.g. has, get) are used to access spectrum parameters, and, in the case where a spectrum generator is connected, the RunToScale function can be called to run the whole spectrum to a different scale. To standardise parameters that should be retrievable from generated spectra for each model, we define a separate SubSpectrumContents object for each model. Each wrapper is then associated with some SubSpectrumContents, and is required to provide all the contents defined therein.

SubSpectrum objects
While the Spectrum interface provides fast and easy access to pole masses, it does not allow access to more detailed spectrum information, such as the values of couplings, mixings, and running mass parameters. For this one must use the SubSpectrum interface. As we covered in Sec. 2.3.2, the simplest way to ob-tain a SubSpectrum object is to extract it from a host Spectrum object using the get_LE()/get_HE() functions or clone_LE()/clone_HE() functions: const SubSpectrum& mssm = spec.get_HE(); unique_ptr<SubSpectrum> copy = spec.clone_HE(); Following the same pattern as the Spectrum interface, spectrum information can be accessed via the SubSpectrum interface using a get function, for which the arguments are a tag, a string identifier, and zero to two indices: using namespace Par; double mass_su1 = mssm.get(Pole_Mass,"∼u",1); Let us now examine the contents of the SubSpectrum interface class. This class is a virtual interface for derived classes individually designed to wrap the data from external spectrum generator codes. The virtualisation allows module writers to treat spectrum data in their own module functions in a generic way, with the selection of the actual spectrum generator to use being a decision left to runtime. The wrapping of external codes is then decoupled from the use of the data required by module functions, meaning that replacement wrappers can be written for new external spectrum generator codes, and the new wrapper activated, without necessitating any modification of downstream module functions.
For reference, a simplified outline of the SubSpectrum class is shown in Figure 6. For the precise class declaration please see the file Elements/include/gambit/Elements/subspectrum.hpp in the GAMBIT source tree.

Getters, checkers and setters
Like the Spectrum interface, the SubSpectrum interface also provides a has function, which accepts the same input arguments as the get function, and which returns a bool that indicates whether or not a given SubSpectrum object contains the requested parameter: bool has_mass_su1 = mssm.has(Pole_Mass,"∼u",1); Unlike the Spectrum interface, the SubSpectrum interface also provides set functions. These can be used to modify or overwrite parameter data contained in the SubSpectrum object. For example, double newmass = 500; // Units generally in double newmass2 = 600; // powers of GeV mssm.set(Pole_Mass,newmass,"∼u",1); mssm.set_override(Pole_Mass,newmass2,"∼u",1); The set and the set_override functions differ in an important way. The set function, if permitted by the underlying wrapper class, will directly change a value in the underlying spectrum generator code, and thus may, for example, impact renormalisation group running. On the other hand, the set_override function does not see the underlying spectrum generator code at all, and will simply store a replacement value that will be preferentially retrieved when using the get function. This can be useful for storing, for example, precision calculations of pole masses in such a way that they can be seamlessly used by other GAMBIT module functions, via the standard interface. We emphasise that the set_override value will take ultimate priority, and in fact the underlying value known to the spectrum generator code will become masked, and irretrievable. So in the example above, a call to mssm.get(Pole_Mass,"∼u",1) after the call to mssm.set_override would retrieve the value 600.
There is also an optional final bool argument to set_override (not shown in Figure 6), which indicates whether the function is permitted to add to the contents of a SubSpectrum. The default is false, leading to an error if a user tries to override a part of the SubSpectrum contents that does not already exist. By setting this true, one can insert arbitrary new information into the SubSpectrum, which will become retrievable with the get method. This functionality should be used sparingly, as the writers of other module functions generally will not be aware of the existence of the additional data. If you find yourself wanting to add large amounts of new data to a spectrum, it is generally better to define a new SpectrumContents, and create a new spectrum wrapper class that conforms to this expanded contents set (see Sec. F.1).

RGE running
The SubSpectrum interface provides the main connection to RGE running facilities in any underlying spectrumgenerator code; the RunBothToScale function of the Spectrum interface is just a thin wrapper over this. These facilities do not exist in all SubSpectrum objects, so one must make sure, via the GAMBIT dependency resolution process, that an appropriate spectrum is supplied to module functions that require running facilities. If a module function attempts to call any of the following functions when no underlying RGEs are connected, then SpecBit will throw an error and terminate.
The two most fundamental RGE-related functions are the RunToScale and GetScale functions. The latter function takes no arguments and simply returns the energy scale at which any running parameters accessible by the get function are defined, for example couplings and running masses. The RunToScale function accepts a new energy scale as input, and when called will run the underlying RGE code and recalculate the running parameters.
The precise behaviour of the RunToScale function can be altered via the optional integer 'behaviour' argument. This controls how attempts to run parameters beyond the known accurate range of the underlying RGEs are handled. There exist two sets of limits in the SubSpectrum interface; upper and lower "soft" limits, and upper and lower "hard" limits. These specify the energy ranges that affect the action of the behaviour switch, and can be checked with the following functions: Note that when writing a new SubSpectrum wrapper these functions need to be overridden to specify valid ranges of RGE evolution, otherwise the default limits (which effectively allow unlimited running) will be used.
The "behaviour" flag acts in conjunction with the hard and soft limits according to the following rules: behaviour = 0 (default) -If running beyond soft limit requested, halt at soft limit (parameter evolution will simply stop silently at this limit). behaviour = 1 -If running beyond soft limit requested, throw warning; beyond hard limit, throw error. behaviour = <anything else> -Ignore limits and attempt running to requested scale (errors may still be thrown by the backend spectrum generator code if it sets its own hard-coded limits).

SLHAea output
Strictly speaking, the SubSpectrum objects do not rely on SLHA, and can interface to the spectrum codes for any model. However, SLHA is an existing standard for the MSSM, so we provide the spectrum in this standard. The SubSpectrum interface thus provides three functions that convert SubSpectrum contents into SLHA2 format, making use of the SLHAea C++ library 6 : void writeSLHAfile(int slha_version, const str& fname) Writes out the contents of the SubSpectrum in SLHA 6 http://fthomas.github.io/slhaea/ format with the name fname. If slha_version is set to 1 then it will conform to SLHA1 standard, while if it is set to 2 then it will conform to the SLHA2 standard. Ideally this should only be used for debugging purposes, as minimising disk access is preferable in high-performance computing applications. However, some backend codes do not have proper APIs and can only be run via the command line, or have otherwise hard-coded requirements for input SLHA files. SLHAstruct getSLHAea(int slha_version) Writes the contents of the SubSpectrum into an SLHAea::Coll object (for which SLHAstruct is an alias). This is an inmemory string-based representation of an SLHA file, which is in the SLHA standard SLHA_version while negating any need to write to disk. Block contents can be easily accessed via [] operators (see the SLHAea documentation for examples). void add_to_SLHAea(SLHAstruct&) Similar to getSLHAea, however the spectrum content is added to a preexisting SLHAstruct, passed in by reference via the input argument and modified in place.
SLHAea is very flexible and can be used to also create SLHA-like structures for other models. Therefore this may also be exploited for non-minimal SUSY models and non-SUSY models too, though other than the MSSM and NMSSM there is no existing standard to fix the precise form and ensure that other codes will understand the results. However, this feature is also rather useful for auto-generated programs, where one can construct the interface between the different generated codes to work with new SLHA blocks.

The SMInputs class
The three data members of a Spectrum object are two SubSpectrum objects, plus an SMInputs object. The latter object does not strictly contain information about the calculated spectrum, and the data it contains is not accessed via the get function of the host Spectrum object. Rather, it contains information that in most cases was used as input to the external spectrum generator code. It is a simple struct that directly mirrors the SMINPUTS, VCKMIN and UPMNSIN blocks defined by the SLHA2 standard (see Figure 7). Because SMInputs is a struct, all data members (the elements of the SLHA2 blocks) have public access, and so parameter values can be accessed directly (see Figure  7 for the variable names). The struct contains only two functions, getSLHAea and add_to_SLHAea. The former simply returns the contents of the struct as an SLHAea object (which can be easily written to disk, for example), while the second adds the contents of the struct to an existing SLHAea object. The latter feature is mainly used internally by Spectrum objects to construct a full SLHA2 representation of their contents.

Extra overloads for get/set/has functions
In the discussions of the interface for Spectrum and SubSpectrum objects, we have mentioned only one set of function signatures for the get, set and has functions used to access spectrum data, however several overloads of these functions are available. Most of these are designed to streamline the interaction between Spectrum objects and the GAMBIT particle database (see Sec 10.1 of Ref. [10]), and they allow particle data to be retrieved according to PDG codes and, where appropriate, antiparticle string names (see Models/src/particle_database.cpp for the full set of definitions known to GAMBIT).
For Spectrum objects, the extra function signatures are as follows: When PDG codes + context integers are used with these functions to specify a particle, the GAMBIT particle database will be used to translate the code into a string identifier, which will then be checked against the string keys defined in the underlying SubSpectrum wrapper to do the data lookup. When writing new wrappers it is therefore important to ensure that particle string names match the strings registered in the GAMBIT particle database.
Similarly, whenever any of these interface functions are called with a string name (potentially also with one index), the underlying wrapper will first be searched for a match to that string. If a match cannot be found, then the GAMBIT particle database will be used to try to translate the string name into a short name plus an index (e.g. "∼chi_0" will be translated to the pair ("∼chi",0)) or vice versa. If there is still no match found with any spectrum contents, the GAMBIT particle database will again be used to try to translate the string name into the string for a matching antiparticle (e.g. "∼chi+_0" will be translated to "∼chi-_0") and the search will be tried again (including translating from string name to short name plus index and vice versa). If no match can be found, an error will be thrown and the scan aborted (except in the case of the has functions, which will simply return false).

Adding support for new models and/or codes
Adding spectrum generators for models not shipped with GAMBIT is an inevitably technical task, as it involves writing an interface to new code. We have however endeavoured to make the procedure as orderly as possible.
To spare more casual users the full technical details, we have placed these in Appendix F. Here we give only a brief overview of the procedure, along with a checklist of tasks to be completed and links to the sections in the Appendix that explain the corresponding details. Some of the jargon used in the checklist has not been necessary to introduce up until now, but rather than explain it here we instead refer readers to the full details in the Appendix.
The bulk of the work is in writing a SubSpectrum wrapper class, which will wrap the new spectrum generator and be usable as we describe in Sec. 2.3. There are also several peripheral tasks to be completed, such as writing a module function to properly construct the new SubSpectrum and package it into a Spectrum object, potentially writing new SubSpectrumContents definitions (which enforce consistency between wrappers that ostensibly share a physics model), and making additions to the GAMBIT model hierarchy (discussed in Ref. [10]). These required tasks are: 1. Choose a SubSpectrumContents definition to which the new wrapper will conform, or write a new one (Sec. F.1). To illustrate the process concretely, we provide a full worked example of the writing of a complicated wrapper for a FlexibleSUSY spectrum generator in Appendix G.

Advanced spectrum usage: vacuum stability
In this Section, we present a sophisticated example of how the Spectrum object and SpecBit module functions can be used, where we calculate stability of the electroweak vacuum for simple Higgs sectors. These functions are included in SpecBit 1.0.0 and will be used in a follow-up to the study of the scalar singlet dark matter global fit [17]. To explain what these module functions do and why, we first present a brief review of vacuum stability. We then give the details of the corresponding observable and likelihood functions implemented in GAMBIT.
The stability of the SM electroweak vacuum in the absence of new physics below the Planck scale (M Pl ≈ 1.2 × 10 19 GeV) has been studied extensively [56][57][58][59][60][61][62][63][64][65][66][67][68] (see also references within, and Refs. [69][70][71]). Although the first bounds on the Higgs mass from vacuum stability were estimated over three decades ago [72,73], the experimentally measured 125 GeV Higgs [74,75] has created increased interest in this topic, as this value results in the existence of a high-energy true vacuum state. If the Higgs field tunnelled from the electroweak vacuum to the global minimum at any point in space, a bubble of low-energy vacuum would form, a process known as bubble nucleation. This bubble would propagate outwards at very nearly the speed of light [56], converting all space in its future light-cone into the true vacuum state, having catastrophic results.
The lifetime of the SM electroweak vacuum does exceed the lifetime of the Universe. It thus exists in a meta-stable state, and the SM is not inconsistent with reality. Yet there is still a non-zero probability of decay, such that a transition to the true vacuum is at some point inevitable.
If new physics does exist below the Planck scale, this may dramatically decrease or increase the lifetime of the false vacuum. For example, as shown in Refs. [62,76,77], using effective couplings suppressed by 1/M Pl , high-energy physics could indeed reduce the lifetime of the electroweak false vacuum down to a fraction of a second. However, even the most minimal extensions to the SM can have a notable impact on the stability of the electroweak vacuum.
Examples of the most recent calculations of the vacuum stability of the SM are those of Ref. [78] and Ref. [68], using three-loop RGEs and two-loop threshold corrections to the quartic Higgs coupling λ at the weak scale. This is referred to as next-to-next-to leading order (NNLO). Before these results the state of the art was next-to leading order (NLO) using two-loop RGEs and one-loop threshold corrections at the weak scale [79][80][81][82][83][84][85][86][87]. These vacuum stability studies were however performed in a flat background spacetime. Studies have also been done with the SM coupled to gravity in a curved spacetime [88,89], motivated by the possibility that the region of the true minimum (typically > 10 10 GeV) may be influenced by the gravitational field.
There have been numerous studies of vacuum stability in BSM theories, such as scalar extensions of the Higgs sector [53,54,[90][91][92][93][94][95]. In these scenarios the addi-tion of a scalar particle can have a positive contribution to the running of the quartic Higgs coupling, which for certain values of the new couplings can completely remove the high energy global minimum [54]. See Ref.
[17] for a discussion of vacuum stability in a scalar singlet-extended Higgs-sector model.
Within GAMBIT we currently use FlexibleSUSY [18] to solve the RGEs up to the Planck scale. For a generic model we are limited to two-loop RGEs and one-loop threshold corrections when using SARAH output. However, higher-order RGEs may be added where these are known, such as in the SM. Therefore even for the scalar singlet model, where we use generated RGEs, this is comparable to the precision of the most recent studies [54,94]. Additionally, as with the GAMBIT spectrum object itself, the vacuum stability module functions may be used in conjunction with any spectrum generator with RGE running functionality that can be interfaced to GAMBIT.
The details of the likelihood function and the relevant physics used to classify a parameter point as either stable, meta-stable or unstable are presented in the following section.

Likelihood details
In this section we outline the details of the derivation leading to the likelihood function used in the GAMBIT vacuum stability module functions. The electroweak vacuum obtains a high energy global minimum when the running quartic Higgs coupling becomes negative at a high scale. For a large renormalisation scale µ the Higgs potential can be approximated with an effective potential of the form The scale at which the quartic coupling becomes negative, and thus the size of the potential barrier between the electroweak vacuum and this global minimum, will affect the likelihood of a quantum mechanical tunnelling to this lower-energy state.
If there exists a high-energy global minimum, then there is a non-zero probability of decay to this state, and thus the electroweak vacuum has a finite lifetime.
There are three possible cases for electroweak vacuum stability, outlined below.
-Stable: If λ(µ) > 0 for all µ < M Pl then the electroweak vacuum is the only minimum of the Higgs potential (up to M Pl ) and is therefore absolutely stable.
-Meta-stable: If there exists a µ 0 < M Pl such that λ(µ 0 ) < 0 and the lifetime of the electroweak vacuum state is longer than the age of the universe. -Unstable: If there exists a µ 0 < M Pl such that λ(µ 0 ) < 0 and the lifetime of the electroweak vacuum state is less than the age of the universe.
Due to the shape of the resulting likelihood function, the distinction between meta-stability and instability is very well defined, with the expected lifetime of the electroweak vacuum being extremely sensitive to changes in the model parameters.
To determine the stability of a model we calculate the rate of quantum mechanical tunnelling through an arbitrary potential barrier, following the derivation in Ref. [96] obtain an analogous field-theoretic solution. In this context, a tunnelling event corresponds to a bubble nucleation in which the Higgs field decays to the lower energy state. We invoke the thin-wall approximation [97], in which the equation of motion of the bubble of radius R, nucleated at time t = 0, is |x| From quantum mechanics, one can derive the rate of tunnelling for a particle through a potential barrier, which can then be generalised to a rate of bubble nucleation per unit volume per unit time of where Γ 0 depends on the size of the past light-cone and B is determined by the shape and size of the potential barrier. Consider a scalar field φ. For bubble nucleation in flat 4-dimensional spacetime we can take φ to be a function of ρ = τ 2 + |x| 2 1/2 only (where τ = it), as there exists an O(4) symmetry. B is given by the action [96] satisfying the Euclidean equation of motion for φ, The boundary conditions of this equation of motion can be combined into the single requirement that lim ρ→∞ φ(ρ) = q 0 . The solution to Eq. 4 is known as the bounce solution, as it describes, in Euclidean time, a solution that is in the unstable vacuum at t → −∞, the global vacuum at t = 0, and bounces back to the starting point at t → ∞. The process outlined above for a simple scalar field can be applied to the Higgs field, with the use of the approximation V (h) ≈ λ(h)h 4 /4. If we perform this calculation at tree level, we can obtain an approximation for the coefficient B. Taking λ(h) = λ, λ < 0, we have a solution [81] to the analogue of Eq. 4 as where R is a dimensional factor associated with the size of the bounce. This can be a relevant quantity such as the height of the barrier, or the change in renormalisation scale between adjacent minima; we shall use the latter. From Eq. 3 we then obtain the action The validity of the approximation V ≈ λh 4 /4, with λ < 0, is not immediately evident, as h = 0 is an unstable maximum, and indeed any value of h > 0 is unstable. However, as stated in Ref. [81], the bounce solution is not a constant field configuration, so requires a nonzero kinetic energy. Because of this, the required bounce solution is suppressed even in the absence of a potential barrier, and thus this is still a valid approximation [98]. Now that we have an expression for B we need to determine the pre-factor. The explicit form of this prefactor, Γ 0 , was first calculated in Ref. [99] by taking account of one-loop quantum corrections. Yet because the behaviour of the exponential in Γ dominates to such a large extent, the value of Γ 0 need only be an approximation, so we follow the analysis of Ref. [56] and determine Γ 0 by a dimensional reasoning.
If we take the Planck constant and the speed of light to be unity ( = c = 1), then the rate of decay per unit time per unit volume, Γ/V has units of [length] −4 . Thus Γ 0 must have units of [length] −4 or [energy] 4 . The characteristic scale relevant in this problem is the width or height of the potential barrier, thus we take Γ 0 ≈ Γ/V is the rate of decay per unit time per unit volume. As we are ultimately interested in the probability of the Universe having decayed in our past light cone, we multiply Γ/V by the volume of the past light cone ∼T 4 U , where T U is the age of the universe. 8 Thus we obtain a rate 7 If the minimum value of the running quartic coupling λmin < 0 is achieved after M Pl , but λ(M Pl ) < 0, we take λmin = λ(M Pl ). 8 This can be computed more rigorously for a standard FLRW cosmology, see Ref. [63], however to the level of detail required here this result is equivalent.
and we obtain the total expected decay rate in the past light cone Here λ(Λ B ) is the minimum value of the quartic Higgs coupling, and we have expressed the age of the Universe as T u ≈ e 140 /M Pl . 9 The arguments used to arrive at Eq. 8 were based only on dimensional analysis. Although this quantity is now dimensionless it cannot be immediately interpreted as a probability, instead it should be interpreted as the expected value for the random variable k, where k is the number of decay events that occurred in the time given (in this case ≈ 10 billion years). 10 To model the probability that the Universe has actually decayed in the given time interval, we use a Poisson distribution, Because we want the likelihood that no decay has occurred in our past light-cone, we calculate the probability that k = 0, which is given by The likelihood given by Eq. 9 is typically either extremely small or exactly one, being extremely sensitive to the value of λ(Λ B ). This results in an almost stepfunction transition from a meta-stable to an unstable Universe when actual Lagrangian-level model parameters are varied. The likelihood in Eq. 9 is difficult to study, due to its double exponential behaviour. In GAMBIT we return the log-likelihood, so this likelihood contains some useful information even within the stable/meta-stable region of the parameter space. In most cases though, this function will either entirely rule out a point or contribute zero to the total log-likelihood. This is indeed useful for ruling out large regions of a parameter space quickly.

Code Description
The vacuum stability calculations in GAMBIT are separated into several module functions, shown in Table  3. As renormalisation group running requires a model spectrum object, these include Spectrum dependencies. The spectrum must be capable of running to different energy scales, and as such cannot be a simple container spectrum.
The module function find_min_lambda, which has capability vacuum_stability, runs the quartic Higgs coupling to determine the scale at which it reaches a turning point. This function produces a struct of type dbl_dbl_bool, which contains the expected lifetime of the Universe (double), the scale of the instability (double), and whether or not the couplings remain perturbative up to the scale returned (bool). If a turning point is not reached before some user-specified high scale (set via YAML option set_high_scale, which defaults to M Pl ), the high scale is returned instead of the scale of the turning point. These calculations can be applied to any model with a Higgs field that does not mix with other scalars.
Derived from the results of find_min_lambda is the physical observable -expected_lifetime: the expected lifetime of the electroweak vacuum in years and the likelihood functions -VS_likelihood: the log-likelihood of vacuum decay not having yet occurred, given the current age of the Universe -check_perturb_min_lambda: the perturbativity of the running couplings at the instability scale. To allow this function to be used as a log-likelihood, it returns 0 for models where the couplings remain perturbative up to the instability scale, and invalidates the model otherwise.
In VS_likelihood, the expected age is compared to the current age of the Universe, which is currently hardcoded to 13.6 × 10 9 years. If λ does not have a turning point before M Pl , then the lifetime returned is 1 × 10 300 years.
The function check_perturb_min_lambda is an efficient way to check that the couplings of a theory remain perturbative up to M Pl or some alternative scale, chosen via the option set_high_scale of function find_min_lambda. By performing this check directly in GAMBIT, we avoid any need to rely on similar checks implemented in some (but not all) spectrum generators. The actual calculation is performed simultaneously with the vacuum stability check, and is thus extracted by check_perturb_min_lambda from the output of the vacuum_stability capability. The result is a flag set true only if all spectrum parameters tagged dimensionless are less than (4π) 1/2 at a set of points evenly distributed in log space between the electroweak scale and the smaller of Λ B and the high scale chosen for find_min_lambda. The default is to test 10 points; this can be reset with the YAML option check_perturb_pts of find_min_lambda.
A simple example of the GAMBIT input YAML commands for using the vacuum stability likelihood, and outputting the Higgs pole mass and the expected age of the universe as observables, is given below.

Higgs couplings
In addition to providing particle spectra to the rest of GAMBIT, SpecBit can also produce Higgs couplings, which it provides as a compact HiggsCouplingsTable. These objects carry all decays of neutral and charged Higgs bosons present in the theory under investigation, their CP eigenvalues, a list of all particles in the theory that Higgses can decay invisibly to, values of the SM-normalised effective couplings of each neutral Higgs to W W , ZZ, tt, bb, cc, ss, τ + τ − , µ + µ − , Zγ and hγ, and additional comparison decays of SM Higgses with the same masses as the neutral BSM Higgses. SpecBit obtains all decays from DecayBit. It makes the CP and invisible state identifications within its own functions. The effective couplings can be obtained in two different ways. The first is to build them directly from couplings provided by dedicated routines in external calculators, as in MSSM_higgs_couplings_FH (see Table 2), which retrieves couplings from FeynHiggs (see Sect. 3.1.3 for details of this interface). Alternatively, they can be estimated on the basis of the partial width approximation, where the squared effective coupling is taken as the ratio of the partial widths to the relevant final state of the BSM and equivalent SM-like Higgs. Examples of the partial width approach are MSSM_higgs_couplings_pwid (Table 2) and SingletDM_higgs_couplings_pwid (Table 3).

DecayBit
All particle decay widths used in GAMBIT are provided by the dedicated module DecayBit. These range from known decays of SM particles, to modifications of SM decay widths and branching fractions due to new physics (e.g. top and Higgs decays), to decay widths of new particles present in BSM theories. These data are needed by ColliderBit [13] for predicting Higgs decay rates and calculating corresponding LHC Higgs likelihoods. They also allow ColliderBit to connect predictions for LHC production of BSM particles with observable final states, by simulating particle collisions and subsequent decays of the particles produced in the initial interaction. The same data are required by DarkBit [12] to simulate cascade decays of final states produced in dark matter annihilation and decay, allowing it to make predictions for cosmic ray fluxes in dark matter indirect detection experiments. Decay widths of SM and BSM particles are also important for predicting the relic density of dark matter, through their participation in thermal freeze-out, or if dark matter is produced non-thermally from decays of some heavy resonance.
The overall structure of DecayBit is outlined in Fig.  8, showing the modular calculation of each particle's decays, and their combination into a complete DecayTable for the model under investigation. Here we begin by describing the decays currently implemented in DecayBit (Sec. 3.1), including the details of the calculations, input data and external codes it employs. We then go through the specific module functions it offers (Sec. 3.2), and the details of the code (Sec. 3.3), covering the associated data structures and utility functions, how to add support for new channels and models, and an explicit worked example of the computation of Higgs boson decays.

Standard Model
DecayBit contains total widths and associated uncertainties for decays of the W , Z, t, b, τ and µ, as well as their antiparticles, and the most common mesons π 0 , π ± , η, ρ 0 , ρ ± and ω. It also includes partial widths to all observed, distinct final states for W , Z, t, b, τ , µ, π 0 and π ± , including many-body final states in many cases (although we do not include such states where the branching fraction would otherwise be contained within a lower-multiplicity channel). We take these data directly from the Particle Data Group compilation of experimental results [51].
We group together hadronic channels for W and Z decays, assigning them a final state consisting of two generic hadrons. Note that 'hadron' is recognised as a generic particle by the GAMBIT particle database; see Sec 10.1 of Ref. [10].
DecayBit also features various calculations of Higgs boson partial and total decay widths. For a pure SM Higgs, the user can choose to either calculate decays at the predicted value of the Higgs mass with FeynHiggs, or extract them from precomputed tables contained in DecayBit. The latter are the predicted SM Higgs partial and total widths given as a function of Higgs mass in Ref. [101]. Using these tables wherever possible is preferable to simply calling HDECAY, as they include additional higher-order corrections for 4-body final states.

Scalar singlet
When m S < m h /2, the Higgs can decay to two singlet scalars. The partial width for this process is Whenever this channel is kinematically open, DecayBit calculates the partial width and adds it to the existing set of SM Higgs decays, rescaling the decay branching fractions (BFs) and total width obtained for a pure SM Higgs.
The decay h → SS would be entirely invisible at the LHC, so would contribute to the invisible width of the Higgs boson. Assuming SM-like couplings, as is the case for the scalar singlet model, 95% confidence level (CL) upper limits on the Higgs invisible width from LHC and Tevatron data are presently at the level of 19% [102]. We implement this as a complete likelihood function in DecayBit, by interpolating the function shown in Fig. 8 of [102].

MSSM
DecayBit will calculate decays of all sparticles and additional Higgs bosons in the MSSM, including branching fractions to SM and SUSY final states. It can also calculate SUSY corrections to decays of the top quark and the SM-like Higgs.
The decays of each Higgs may come (independently) from backend functions in either HDECAY via SUSY-HIT, or FeynHiggs. For the calculation of effective Higgs couplings by SpecBit (see Sec. 2.6), DecayBit will also compute decays of an SM Higgs with the same mass as any of the MSSM Higgs bosons, using either FeynHiggs or by interpolating its own internal tables. For masses between 80 GeV and 1 TeV, these tables come directly from Ref. [101]. Between 1 GeV and 80 GeV, and between 1 TeV and 16 TeV, we supplement them with additional results that we obtained by running HDECAY 6.51 [26]. For Higgs masses above 1 TeV, in order to obtain finite results we disabled electroweak corrections and decays to ss, bb, gg and µ + µ − (which are negligible at large masses, where W + W − , ZZ and tt final states dominate).
Top quark decays are available only by calling Feyn-Higgs. Sparticle decays are available only from SDECAY via SUSY-HIT. DecayBit can calculate Higgs and top widths for arbitrary MSSM-63 models when employing FeynHiggs, but is limited to MSSM-20 models for sparticle and Higgs decays when using SUSY-HIT, due to the fact that SUSY-HIT is not SLHA2 compliant.
We backend SUSY-HIT 1.5 in a similar fashion to other codes, but we patch it via the GAMBIT build system when it is automatically downloaded, in order to make a number of crucial modifications. Here we go through the most important of these.
SUSY-HIT is written as a standalone program, whereas GAMBIT requires callable functions in a shared library. We therefore convert the main program of SUSY-HIT to a function, and modify the makefile to produce a shared library.
Interfacing to HDECAY and SDECAY via SUSY-HIT means that we retain the standard HDECAY and SDECAY options chosen there: multi-body decays, loopinduced decays and 1-loop QCD corrections to 2-body decays are all enabled, and any running masses and couplings are calculated at the EWSB scale in the DR scheme. However, we directly inject the bottom quark pole mass calculated by SpecBit into SUSY-HIT, rather than relying on its internal pole mass calculation.
For some rather specific models, the 1-loop QCD corrections ∆ 1 to the widths of the lighter sparticles are negative and larger than the tree-level results, causing SUSY-HIT to return negative decay widths. We have so far found models where this is the case for gluinos, sbottoms, stops, sups, sdowns, neutralinos and charginos. To rectify this, we have patched SUSY-HIT to implement a correction factor designed to approximately mimic resummation. When the corrections are negative and larger than the tree-level result, instead of obtaining the 1-loop corrected widths as usual in SUSY-HIT we apply the correction as However, this no longer follows when Γ tot is corrected by Eq. 12. In this case, we simply revert to determining branching fractions from the tree-level partial and total widths, but retain the more accurate total width obtained by Eq. 12 as the total decay width to be used elsewhere in GAMBIT.
Correcting negative widths with Eq. 12 has a downside: the total width of a sparticle becomes discontinuous in the model parameters at ∆ 1 = −1. Whilst this is still an obvious improvement over having non-physical results for ∆ 1 < −1, it is not ideal. An alternative treatment would be to use Eq. 12 whenever ∆ 1 < 0, which would ensure continuity, and to apply the correction to each partial width individually, allowing one-loop branching fractions to still be calculated when ∆ 1 < 0. However, this would modify the standard behaviour of SUSY-HIT in the regime −1 < ∆ 1 < 0. In the interests of keeping our modifications to external codes simple and minimal, we therefore consider it preferable to pay

Reference_SM_Higgs_decay_rates lnL_Higgs_invWidth lnL_Higgs_invWidth_SMlike(double):
Computes the log-likelihood of the Higgs invisible width. Table 5: Scalar singlet functions provided by DecayBit. Neither of these functions has any backend requirements nor options.

Higgs_decay_rates
the price of a discontinuous width, restricting our modifications to changes that only have an impact when SUSY-HIT would otherwise return unphysical results.
A related issue is that at high SUSY-breaking scales, the M S correction to the tree-level t and b pole masses required for calculating Higgs couplings in SUSY-HIT also becomes large, leading to numerical instabilities. To remedy this, we import the improved running M S expressions for m t and m b from HDECAY 6.51 to SUSY-HIT.
Finally, we also implemented a diskless generalisation of the SLHA interface to SUSY-HIT. This avoids the need to read and write files to disk, allowing GAMBIT spectrum contents to be directly injected into SUSY-HIT common blocks in SLHA format, and decay data to be extracted directly from other common blocks and easily converted to the native GAMBIT DecayTable format (described in Sec. 3.3).
We connect to and use FeynHiggs as a regular GAMBIT backend, without significant 11 modification. GAMBIT supports versions 2.11.2 and later. We adopt the settings for the real MSSM recommended in the FeynHiggs documentation, but restrict neutral Higgs mixing to the CP -even states, for compatibility with HiggsBounds/HiggsSignals. For version 2.12.0 we set loglevel = 3 in order to include NNLL corrections.

SUSY-HIT
snuX_taul_decay_rates snuX_taul_decays(DecayTable::Entry): Computes decays of theντ L sneutrinos.  We pass MSSM models to FeynHiggs using the function FHSetPara, directly setting the SM parameters in SLHA2 parameterisation, the CP -odd Higgs pole mass, the µ parameter at the SUSY scale and the DR soft terms (sfermion mass parameter matrices, gaugino mass parameters and trilinear couplings, also at the SUSY scale), from the outputs of the Spectrum objects provided by SpecBit.

Available functions and options
Here we give a compact census of the module functions in DecayBit. These are also laid out in Tables 4-7.
Antiparticle equivalents of all functions described below are automatically created from the decays of their particle partners, assuming CP symmetry. This need not be assumed when adding additional states or channels to DecayBit.

Standard Model
The SM particle decay functions are listed in Table 4 Table 2). The result of one or the other of these two functions is presented to the rest of GAMBIT as the definitive Higgs_decay_rates, by function SM_Higgs_decays.
The two functions are also used by other module functions to compare BSM Higgs decays to equivalent SM decays.

Scalar singlet
We highlight the parts of DecayBit specific to the scalar singlet model in Table 5. Only the Higgs_decay_rates are modified in this model. The function SingletDM_Higgs_decays simply piggybacks off the SM Higgs decay calculation, and then rescales the total width and branching fractions to accommodate the additional h → SS contribution. The SM-like Higgs invisible width likelihood is provided by the function lnL_Higgs_invWidth_SMlike, which has capability lnL_Higgs_invWidth and depends on the Higgs_decay_rates. Table 6 lists the full set of MSSM decay functions contained in DecayBit. The naming scheme follows a similar structure of the SM decays, with function names particle_designation_decays and capabilities particle_designation_decay_rates.

MSSM
Here particle_designation includes the particle name and any relevant indices: plus or minus for charged bosons, bar for sfermionic antiparticles, mass ordering indices for bosons, neutralinos and charginos and third-generation sfermions, and weak eigenstate indicators l or r for first and second-generation sfermions. Despite the fact that GAMBIT works exclusively in SLHA2 mass-ordered basis internally, we designate first and second-generation sfermion decays in the weak eigenstate basis because SUSY-HIT is not SLHA2-compliant, and left-right mixing is typically small in the first two generations in the MSSM anyway. This will probably be generalised somewhat further in future revisions of DecayBit. We discuss the conversion between mass, weak and family-mass eigenstates in more detail in Sec. 3.3.2.
DecayBit can compute MSSM Higgs sector decays using either SUSY-HIT or FeynHiggs; it includes separate module functions for each of these options. It also includes additional MSSM-mediated top decays, calculable via FeynHiggs. Similar to the reference functions for the SM-like Higgs described in Sec. 3.2.1, De-cayBit also provides functions to compute reference decays of an SM Higgs with the same mass as the A 0 (capability Reference_SM_A0_decay_rates) and the 'other', non SM-like, CP -even neutral MSSM Higgs (usually -but not always -the heavy Higgs; capability Reference_SM_other_Higgs_decay_rates). These are neccessary for constructing effective Higgs couplings from partial widths with SpecBit (see Sec. 2.6).

Collectors and helpers
DecayBit also ships with a number of other useful module functions, given in Table 7. The most important are the collector functions all_decays and all_decays_from_SLHA.
These both have capability decay_rates and return a complete GAMBIT DecayTable object (see Sec. 3.3.1 below). This is a single object with all decays of all known particles in the models under investigation, to be used by the rest of GAMBIT in subsequent physics calculations. The first of these functions constructs the DecayTable from the results of the relevant individual decay module functions in DecayBit. The second is a convenience/debug tool that simply bypasses the rest of DecayBit, reading some pre-computed decays from a user-supplied SLHA file and presenting them to the rest of the code as if they had been computed by DecayBit.
The other two module functions have to do with the conversion between sfermion weak/family and mass eigenstates. The function check_first_sec_gen_mixing checks whether first or second generation weak eignestates are mixed beyond the level specified by the option gauge_mixing_tolerance, returning zero if all states are sufficiently pure and an integer between 1 and 6 otherwise, depending on which eigenstate shows the mixing. This function is mostly intended for saving information about mixing for later inspection. The actual sparticle decay functions all have dependencies on the SLHA_pseudonyms, which list the sparticle mass eigenstates containing the largest admixture of each weak eigenstate; the function that returns this list is get_mass_es_pseudonyms.
The underlying weak/family-mass eigenstate conversion routines are discussed in more detail in Sec. 3.3.2.

The DecayTable
GAMBIT includes a specific DecayTable container class for storing and passing decay information to the rest of the code. A DecayTable contains one or more instances of the subclass DecayTable::Entry, each of which holds the full decay information for a single particle. The primary job of DecayBit is to construct a DecayTable::Entry for each particle in a given model, and then to put them together into a complete DecayTable that other parts of the code can access.
Each DecayTable::Entry includes the width of the particle, the positive and negative uncertainties on that width, the branching fractions to different decay channels, the uncertainties on the branching fractions, and information about which code(s) (GAMBIT, a backend or some combination) was responsible for computing the decay information. The width and branching fraction uncertainties are not currently used in any calculations in GAMBIT, and are only set for SM particles (where they are known experimentally from the Particle Data Group [51]), but they may be set more completely and used in more advanced likelihood functions in the future, for example in Higgs physics. Widths can be set and accessed directly with the width_in_GeV field. Decay channels and branching fractions are added and retrieved with set_BF() and BF(). The following simple example for the decays of the W + boson shows these in action: // Declare a set of decays for the // W + DecayTable::Entry Wplus; // Set the total width Wplus.width_in_GeV = 2.085; // Set branching fraction and uncertainty // for W → e νe Wplus.set_BF(0.1071, 0.0016, "e+", "nu_e"); // Find the partial width for W → e νe double pw = Wplus.width_in_GeV * Wplus.BF("e+","nu_e"); The functions set_BF() and BF() each have five different overloads, allowing the decay channel to be specified with any of the particle identifiers recognised by the GAMBIT particle database (see Sec. 10.1 in Ref. [10]). These are: 1. a simple argument list of strings corresponding to the final state particles' long names (as in the example above) 2. an argument list of alternating PDG codes and context integers, 3. an argument list of alternating short name strings and integer particle indices, 4. a std::vector of long name strings, and 5. a std::vector<std::pair<int,int>> of PDG codecontext integer pairs.
All of these forms allow an arbitrary number of final state particles to be specified for a given decay channel. The DecayTable::Entry subclass also has a simple convenience method has_channel() that takes channels specified in any of the above forms, and returns a bool indicating whether the Entry already contains a field corresponding to that decay channel or not. The individual entries in a full DecayTable can be accessed using semantics reminscent of a regular C++ std::map, where the key is the decaying particle, referred to by any of the three possible designators recognised by the GAMBIT particle database. For example, // Declare a DecayTable DecayTable BosonDecays // Add entries (here are some I prepared earlier) BosonDecays("W+") = Wplus; BosonDecays("W-") = Wminus; BosonDecays("Z0") = Z; BosonDecays("h0_1") = SM_higgs; // Get the branching fraction for W → e νe double enueBF = BosonDecays("W+").BF("e+","nu_e"); // Get underlying Entry for Higgs // six different ways DecayTable::Entry h; h = BosonDecays("h0_1"); h = BosonDecays(25, 0); h = BosonDecays("h0", 1); h = BosonDecays.at("h0_1"); h = BosonDecays.at(25, 0); h = BosonDecays.at("h0", 1); In this example, a DecayTable is created for holding the decays of SM bosons. The decays of the individual channels, each an instance of the DecayTable::Entry subclass, are then added to the table with the () method. A branching fraction for W → e ν e is then extracted directly from the DecayTable. The entire entry for Higgs decays is then extracted in six different ways, showing the standard bracket methods, the constant at access method (reminiscent of regular C++ STL container access), and the three different ways of referring to the particle whose decays are desired from the table.
The DecayTable machinery integrates seamlessly with SLHA and its generalisations, via SLHAstruct SLHAea objects. Individual DecayTable::Entry instances can be constructed directly from SLHA DECAY blocks provided in SLHAea::Block format, and made to emit them with the method getSLHAea_block. Likewise, a complete DecayTable can be constructed directly from an SLHA file or an SLHAstruct object containing DECAY blocks, and can be output to either format with writeSLHA and getSLHAea.
Full documentation of the programming interface provided by the DecayTable and DecayTable::Entry classes can be found in Elements/include/gambit/Elements/decaytable.hpp.

Utilities
The module functions check_first_sec_gen_mixing and get_mass_es_pseudonyms (Table 7) rely on a series of helper functions that probe the sfermionic MSSM mixing matrices provided by SpecBit, in order to determine the gauge composition of a given mass or family-specific mass eigenstate, or the converse for gauge eigenstates. These functions are most relevant for DecayBit and backends such as SUSY-HIT, but are needed also by Col-liderBit and DarkBit. They therefore ship as part of the main GAMBIT distribution, in the Elements directory. Functions can be found here for returning: the mass eigenstate with largest contribution from a given gauge eigenstate, plus full mixing and composition information if required (slhahelp::mass_es_from_gauge_es()) -the gauge eigenstate with largest contribution from a given mass eigenstate, plus full mixing and composition information if required (slhahelp::gague_es_from_mass_es()) -the (6-flavour) mass eigenstate that best matches a given family (2-flavour) mass state, plus full mixing and composition information if required (slhahelp::mass_es_closest_to_family()) -the family (2-flavour) mass state that best matches a given (6-flavour) mass eigenstate, plus full mixing and composition information if required (slhahelp::family_state_closest_to_mass_es()) -the two (6-flavour) mass eigenstates that best match a requested generation, as well as the implied (2flavour) family mixing matrix between the family states. (slhahelp::family_state_mix_matrix())

Full documentation can be found in
Elements/include/gambit/Elements/mssm_slhahelp.hpp.
The DecayBit module proper contains just two utilities of note that are not actual module functions (in the GAMBIT sense).
The first is the class DecayBit::mass_es_pseudonyms, which is used as the return type for get_mass_es_pseudonyms after obtaining the necessary content from the slhahelp functions. This is essentially just a container of named strings, one for each SLHA1 MSSM weak/family eigenstate sfermion.
The second is the function DecayBit::CP_conjugate(), which takes as its lone argument a DecayTable::Entry and returns a corresponding DecayTable::Entry for the CP conjugate of the original particle, assuming CP symmetry. The X=bar and Y=minus variants of the module functions in Tables 4 and 6 consist of nothing more than a call to this function, using their dependencies on the corresponding X=empty and Y=plus variants.

A worked example: SM-like Higgs decays
To illustrate some of the workings of DecayBit more concretely, here we go through the example of decays of the Higgs in detail.
The rollcall header DecayBit_rollcall.hpp declares a series of functions that are each capable of generating a DecayTable::Entry containing the width and branching fractions for Higgs decays: The first of the these, Ref_SM_Higgs_decays_table, calculates the width and decay branching fractions of the pure SM Higgs, interpolating in the tables of Ref. [101] using the Higgs pole mass (provided by another module function as a dependency). This decay information is presented to the rest of the code with capability Reference_SM_Higgs_decay_rates, to allow other functions to declare that they depend on knowing the pure SM decays as reference quantities, even if the actual decays of the Higgs in the model being scanned turn out to be different to the SM case. The function SM_Higgs_decays has a dependency on this reference set of decay information, and simply presents it 'as is' to the rest of the code as the definitive decays of the Higgs, as indicated by the capability Higgs_decay_rates. The module function SingletDM_Higgs_decays gives the Higgs decays in the case of an SM-like Higgs that can also decay to two singlets. This function also piggybacks off the SM reference decay information, simply adding in the h → SS decay and rescaling the SM branching fractions. MSSM_h0_1_decays returns decays of the lightest MSSM Higgs.
The latter two functions have corresponding ALLOW_MODELS declarations, whereas SM_Higgs_decays is considered to be universally compatible with all models, and will be used as a fallback in any case where a modelspecific Higgs decay function does not exist, i.e. assuming that the Higgs is unaltered from the SM. The singlet variant also has a dependency on the SingletDM_spectrum, which it uses to extract the S mass and coupling in order to compute the invisible width.
The MSSM variant has no dependency on a Spectrum object, because it instead relies on BACKEND_REQuirements, which it demands to be filled from its only declared BACKEND_OPTION, SUSY-HIT. A spectrum object is not needed because the backend initialisation function for SUSY-HIT has this dependency itself, as that is where the pole masses, couplings and mixings are actually fed into SUSY-HIT. The function MSSM_h0_1_decays does have a dependency on the SLHA_pseudonyms however, as these are needed to interpret the weak/family eigenstate decay widths returned by SUSY-HIT, and to place them into the resulting DecayTable::Entry in terms of mass eigenstates.
The macro calls illustrated in these declarations are described in more detail in the full GAMBIT paper [10].
The SM reference function computes the total width and branching fractions of the Higgs as follows: std::stringstream msg; msg << "Requested Higgs virtuality is " << mh << "; allowed range is " << minmass << "--" << maxmass << " GeV."; invalid_point().raise(msg.str()); } // Set the contents of the Entry result.calculator = "GAMBIT::DecayBit"; result.calculator_version = gambit_version; result.width_in_GeV = virtual_SMHiggs_widths ("Gamma",mh); result.set_BF(virtual_SMHiggs_widths("bb",mh), 0.0, "b", "bbar"); ... } First the Higgs mass is extracted from the mh dependency, and checked against the allowed minimum and maximum Higgs mass. These bounds are used to ensure that the SM-like Higgs is not outside the range over which the values covered by the tables of Ref. [101] are most reliable. The calculator and version are then set, and the total width is extracted from the Ref. [101] tables with the function virtual_SMHiggs_widths (found in Elements/src/virtual_higgs.cpp). Finally, branching fractions to different final states are set by calls to set_BF(), in a similar fashion to the first example of Sec. 3.3.1.
The body of SM_Higgs_decays simply rebrands this information from Reference_SM_Higgs_decay_rates to Higgs_decay_rates, indicating to the rest of GAMBIT that it indeed represents the true Higgs decays in this case, and should be used as such in later calculations: Here we see the generic dependency upon Higgs_decay_rates; regardless of the model being scanned, all_decays requires a DecayTable::Entry that describes the decays of the Higgs boson. Depending on which model is actually being scanned, this can end up being fulfilled by any one of the three functions with capability Higgs_decay_rates described above. The function all_decays also has a MODEL_CONDITIONAL_DEPENDENCY on the set of decays of each BSM particle, such as the second CP -even Higgs in the MSSM. The function body of all_decays simply takes all the entries that it depends on, and adds them to a single DecayTable object, along the lines of the second example in Sec. 3.3.1.

Adding support for new models and programs
Supporting a new model in DecayBit is a matter of: 1. writing additional module functions that can generate DecayTable::Entry objects for the novel field content, 2. adding new module functions that can compute new decays of existing particles in the new model, and Model parameters are used by SpecBit to provide masses and couplings to PrecisionBit, and also by PrecisionBit directly, to compute most SM nuisance likelihoods. PrecisionBit uses the precision Higgs masses from SpecBit, along with its own precision calculation of the W mass, to convert the spectrum provided by SpecBit into a precision-updated version. It then uses the improved spectrum for calculating likelihoods associated with the Higgs and W masses, as well as a suite of electroweak precision observables (EWPO). It then feeds the EWPO predictions to its own experiment-based EWPO likelihoods.

adding dependencies on these new functions to
all_decays.
Adding a new backend capable of calculating decay rates and branching fractions is essentially the same as adding any new backend, and should follow the recipe presented in the main GAMBIT paper [10]. After this, using the new backend from DecayBit requires writing new module decay functions designed to take advantage of it, or modifying existing ones to do so. In some rare cases, the backended functions might have the same signature as those of an existing backend, in which case the backend simply needs to be selected at the YAML file level, using the Rules section (see Sec. 6 of [10]).

PrecisionBit
The PrecisionBit module serves a number of related purposes. It provides so-called nuisance likelihoods, which describe the uncertainties on known quantities that are important inputs to other calculations, which one might like to vary within their experimentally-allowed ranges in the course of a scan. This functionality is discussed in Sec. 4.1. Examples include the mass of the top quark and the strong coupling; similar quantities relevant only to dark matter observables are dealt with in DarkBit (e.g. the strange quark content of the proton and the local density of dark matter).
PrecisionBit is also responsible for calculating BSM corrections to precision SM observables, such as the mass of the W boson and the weak mixing angle θ W , and providing correspondingly corrected pole masses and couplings to the rest of GAMBIT wherever appropriate. It also provides likelihood functions for quantifying the agreement of the predicted corrections to the precision observables with experimental data. We discuss Preci-sionBit's ability to compute MSSM precision corrections and likelihoods in Sec. 4.2. Fig. 9 shows how these different components of Pre-cisionBit fit together, and how the two 'halves' of the module are connected via SpecBit and the input SM parameters.

Standard Model nuisances
PrecisionBit provides a number of SM nuisance likelihoods, related to SM couplings, quark masses and the masses of the Z and W bosons. These likelihoods are generally approximated as χ 2 -like functions, comparing experimentally determined central values (and errors) with the SM input values provided by GAMBIT. Whenever they are available and relevant, we add theoretical errors to the experimental ones in quadrature.

Couplings
Likelihoods associated with the Fermi coupling constant (G F ) and fine-structure constant (α em ) are constructed by comparing model input values with those determined from electroweak global fits, with G F = (1.1663787 ± 0.0000006) × 10 −5 GeV −2 and α em (m Z ) −1 = 127.940 ± 0.014 (M S scheme) [51]. Similarly, the strong coupling (α s ) likelihood is calculated using a χ 2 -averaged value of α s (m Z ) = 0.1185 ± 0.0005 (M S scheme) taken from the sub-field of lattice determinations [51]. Uncertainties are taken to correspond to 1 σ confidence intervals with no additional theoretical uncertainties. The PrecisionBit functions providing these likelihoods are summarised in Table 8.

SMINPUTS lnL_alpha_s lnL_alpha_s_chi2(double):
Computes the log-likelihood of the strong coupling constant, αs(m Z ), in the M S scheme.

SMINPUTS lnL_GF lnL_GF_chi2(double):
Computes the log-likelihood of the Fermi coupling constant, G F .

SMINPUTS lnL_light_quark_masses lnL_light_quark_masses_chi2(double):
Computes the joint log-likelihood of the M S masses of the u, d and s quarks at µ = 2 GeV.

SMINPUTS lnL_t_mass lnL_t_mass_chi2(double):
Computes the log-likelihood of the top quark pole mass.

SMINPUTS lnL_Z_mass lnL_Z_mass_chi2(double):
Computes the log-likelihood of the Z boson pole mass

SMINPUTS lnL_W_mass lnL_W_mass_chi2(double):
Computes the log-likelihood of the W boson pole mass.

mw lnL_h_mass lnL_h_mass_chi2(double):
Computes the log-likelihood of the h boson pole mass.
mh Table 8: PrecisionBit module functions providing log-likelihood calculations related to SM couplings and particle masses. All functions depend on SMINPUTS with the exception of the likelihoods of the W and h masses, which depend directly on the calculated W and h masses, respectively. These are available from both SM and BSM spectra, as summarised in Table 9. Except in very simple scenarios, the more detailed Higgs likelihoods of ColliderBit [13] are generally to be preferred to lnL_h_mass_chi2, as they simultaneously take account of correlated constraints on the mass and couplings.
GAMBIT with the corresponding PDG values and uncertainties [51], approximating the likelihoods as χ 2 functions. Similarly, the likelihood corresponding to the top quark mass is evaluated using combined Tevatron and LHC measurements [106]. For the light quark masses a single, joint likelihood is calculated with terms corresponding to the ratio of u and d masses, the ratio of the s mass to the sum of u and d masses, and m s , respectively. Users can optionally provide their own values for these three quantities.
The W and Z boson mass likelihoods are based on mass measurements and corresponding uncertainties from the Tevatron and LEP experiments [51]. These strictly correspond to the mass parameter in a Breit-Wigner function rather than pole masses, but the difference is negligible for almost all purposes. DecayBit also provides a simple Higgs mass likelihood, based on the PDG [51] combination of ATLAS and CMS measurements at the LHC. The module functions providing the W , Z and h mass likelihoods are summarised in Table 8. Unlike the other SM particle masses, the modeldependent W and h mass dependencies are satisfied by functions specific to different model spectrum types, listed in Table 9. In this case, if theoretical uncertainties are available, we add them in quadrature to the experimental uncertainties when evaluating the χ 2 approximation to the mass likelihoods. The more detailed Higgs likelihoods of ColliderBit [13] will usually be more useful than the simple Higgs mass likelihood provided in DecayBit. The latter can offer a lightweight alternative when varying the mass of the SM Higgs as a nuisance parameter in simple scans, as in e.g. Ref. [17].

MSSM precision observables
PrecisionBit allows for users to calculate a collection of precision observables using an array of external code packages. These observables can be accessed directly through module functions, passed as input to likelihood calculations, or used to update particle spectra for use throughout GAMBIT.

External code interfaces
Currently, PrecisionBit makes use of the external code packages FeynHiggs [21], SuperIso [28] and GM2Calc [29] for calculating precision observables in the MSSM. The GAMBIT interface to FeynHiggs was described already in Sec. 3.1.3, and the interface to SuperIso is covered in detail in the FlavBit paper [14].

Capability Function (Return Type): Brief Description
Dependencies mw mw_from_SM_spectrum(triplet<double>): Provides the W mass and its uncertainties from an SM spectrum.

SM_spectrum mw_from_SS_spectrum(triplet<double>):
Provides the W mass and its uncertainties from a scalar singlet model spectrum.

SingletDM_spectrum mw_from_MSSM_spectrum(triplet<double>):
Provides the W mass and its uncertainties from an MSSM spectrum.

MSSM_spectrum mh mh_from_SM_spectrum(triplet<double>):
Provides the h mass and its uncertainties from an SM spectrum.

SM_spectrum mh_from_SS_spectrum(triplet<double>):
Provides the h mass and its uncertainties from a scalar singlet model spectrum.

SingletDM_spectrum mh_from_MSSM_spectrum(triplet<double>):
Provides the mass and associated uncertainties of the most SM-like Higgs boson in the MSSM. The module function GM2C_SUSY calls library routines from GM2Calc to calculate the anomalous magnetic moment of the muon and translates the result and the errors from a µ to (g − 2) µ . GM2Calc can accept parameters and masses given in SLHA conventions. The module function exploits GAMBIT helper functions, discussed in Section 3.3.2, to convert between SLHA and SLHA2, which GAMBIT uses internally to store the spectrum. The inputs from the spectrum object, including the SMINPUTS parameters, are then passed on to the GM2Calc object by the module function GM2C_SUSY. One important subtlety is that GM2Calc also has two non-SLHA inputs, which provide input values for the fine structure constant, α, in the Thompson limit and another α that includes SM fermion contributions to the on-shell photon vacuum polarization. These definitions differ from the M S α(m Z ), the inverse of which appears in entry 1 of the SLHA input block SMINPUTS. By default, GAMBIT uses the default values recommended by the GM2Calc authors. It is possible to overwrite these values using the YAML options GM2Calc_extra_alpha_e_MZ and GM2Calc_extra_alpha_e_thompson_limit.

MSSM_spectrum SMlike_Higgs_PDG_code
PrecisionBit also contains an experimental interface to the proprietary code SUSY-POPE [107,108], but this is not recommended for general use.
The various interface functions to different external code packages within PrecisionBit are summarised in Table 10, and include calculations of backend-specific data-types for input to other functions, accessor functions for specific observables, and functions for updating particle spectra to include precision corrections to the W and Higgs-sector masses.

Electroweak precision observable likelihoods
PrecisionBit includes a collection of functions for using the precision observables calculated by external code packages to calculate likelihoods. Currently, likelihood functions are provided for the effective leptonic weak mixing angle sin 2 θ W,eff , the departure from 1 of the ratio of the Fermi constants implied by neutral and charged currents ∆ρ, and the muon anomalous magnetic moment a µ (also referred to as g − 2 ≡ 2a µ ). For each of these observables, PrecisionBit calculates likelihoods by combining experimental and theoretical uncertainties in quadrature, and constructing likelihood functions based on the difference between the calculated and experimentally-measured values. The module functions that provide these capabilities are listed in Table 11.
For the g − 2 likelihood, we combine errors from the predicted SM contribution in quadrature with the theoretical error arising from the BSM contribution. In the likelihood function, we compare the sum of the SM and BSM contributions to the experimental measurement a µ = (11659208.9 ± 6.3) × 10 −10 [110,111], where the error is the sum in quadrature of statistical (5.4 × 10 −10 ) and systematic (3.3 × 10 −10 ) uncertainties.
Theoretical calculations of the MSSM contribution to a µ can be taken from either GM2Calc, SuperIso or FeynHiggs. While SuperIso and FeynHiggs include essentially the same one and two-loop corrections, and give more or less consistent results, GM2Calc includes additional two-loop contributions and an improved on-shell calculation of the one-loop effects, leading to significant improvements in precision for some parameter combinations. GM2Calc calculates an uncertainty on

SP_PrecisionObs SP_PrecisionObs(double):
Calculate precision observables in the MSSM with SUSY-POPE (experimental). Table 10: PrecisionBit module function interfaces to external code packages for calculation of precision observables. These functions include calculations of backend-specific data-types for input to other functions, retrieval functions for specific observables, and functions for updating particle spectra for use throughout GAMBIT. Computes the SM contribution to g − 2, based on e + e − data.

gm2_SM_tautau(triplet<double>):
Computes the SM contribution to g − 2, based on τ + τ − data. its prediction, using the magnitudes of the two-loop Barr-Zee corrections [112] to estimate the magnitude of neglected higher-order contributions [29]; PrecisionBit adopts this estimate directly when employing GM2Calc.
For FeynHiggs and SuperIso, no such estimate is available. Following the discussion in Ref. [113], we assign a theoretical uncertainty of either 30% or 6 × 10 −10 (whichever is greater) to the values of g −2 that PrecisionBit obtains from SuperIso and FeynHiggs. PrecisionBit's theoretical predictions of sin 2 θ W,eff and ∆ρ in the MSSM are currently provided exclusively by FeynHiggs. We assume a theoretical error of 12 × 10 −5 on sin 2 θ W,eff [114]. The MSSM corrections to both sin 2 θ W,eff and m W are approximately proportional to ∆ρ at tree level. Combined with the fact that ∆ρ 1 in general, to a good approximation the theoretical error on its value can be estimated from the fractional uncertainty on sin 2 θ W,eff or m W ; we therefore adopt a theory error on ∆ρ of where σ theo m W = 10 MeV [114].

Precision-updated MSSM spectrum
As described in Table 10, PrecisionBit includes a number of module function for updating existing MSSM spectra with precision calculations of particle masses: These functions take an unimproved_MSSM_spectrum as a dependency, update its mass spectrum with precision values, and return the improved version as an MSSM_spectrum. The H_W version updates both the Higgssector and W masses, the W version does only the W mass, and make_MSSM_precision_spectrum_none simply transmutes an unimproved_MSSM_spectrum into an MSSM_spectrum without further modification. The precision W and Higgs-sector masses are typically provided via dependencies fulfilled by functions that call Feyn-Higgs (though they can of course in principle come from anywhere). In addition to updating the masses, the make_MSSM_precision_spectrum_H_W() and make_MSSM_prec-ision_spectrum_W() functions also update the associated theoretical uncertainties provided by the original spectrum object. The uncertainty on the pole mass of the W is adopted from its dependency on the precision value of m W . When calculating this with FeynHiggs, we assume an uncertainty of 10 MeV [114] (as discussed in the context of the theory error on ∆ρ in Sec. 4.2.2). 13 For the MSSM Higgs bosons (A 0 , h 0 , H 0 and H ± ), the methods by which PrecisionBit arrives at the central value and uncertainties on the mass are configurable by the user, using the options Higgs_predictions_source and Higgs_predictions_error_method, respectively. For the Higgs_predictions_source, there are three possible values: 1: Update each of the Higgs masses to the value derived from the precision calculator (canonically FeynHiggs). This is the default. Five options are available for assigning theoretical uncertainties to the MSSM Higgs masses. These are based on different combinations of three quantities with unique values for each Higgs state: the error associated with the mass from the original input spectrum (∆ s ), the error determined by the precision calculator (∆ p ), and the difference of central values from the two calculations (∆ g ). All but one of the options for using these quantities assigns them as additional errors beyond a so-called 'range around chosen central' (RACC) calculation. The RACC is defined such that the upper error σ + is the distance from the central value chosen with Higgs_predictions_error_method to the greater of the values predicted by the original spectrum and the precision calculator. The lower error σ − is defined analogously as the distance from the central value chosen to the lesser of the two predictions. With these definitions, the available options for the Higgs_predictions_error_method are: 1: Upper and lower uncertainties both set equal to the sum in quadrature of all three uncertainties: RACC, with ∆ s added to the error associated with the distance to the spectrum generator value (the spectrum-generator 'edge'), and ∆ p added to the precision-calculator edge. This is the default. 3: RACC, with ∆ g /2 added to both σ + and σ − . 4: RACC, with ∆ g /2 added at the spectrum-generator edge, and ∆ p added at the precision-calculator edge. 5: RACC, with ∆ g /2 added at the precision-calculator edge, and ∆ s added at the spectrum-generator edge.
The resulting precision-improved spectrum is provided to the rest of GAMBIT for physics calculations, and in the case of the Higgs bosons, computation of Higgs-mass likelihoods by ColliderBit.

Example YAML files
In yaml_files, we give some example input YAML files that run various functions from SpecBit, DecayBit and PrecisionBit within GAMBIT. The file SpecBit_MSSM.yaml runs a number of tests on the various MSSM spectra, using either FlexibleSUSY or SPheno from within SpecBit. The file SpecBit_vacuum_stability.yaml does a short (10 minutes runtime on a single core) scan of the Higgs and top mass, in order to map out the regions of SM vacuum (meta-)stability and perturbativity. The DecayBit example DecayBit_MSSM20.yaml computes decay widths and branching fractions for a single example point in the MSSM20atQ model [10]. DecayBit_SingletDM does the same for ten randomlyselected parameter combinations in the scalar singlet dark matter model. PrecisionBit_MSSM20.yaml computes all precision observables and likelihoods contained in PrecisionBit, for the same example MSSM20atQ point as DecayBit_MSSM20.yaml.

3-BIT-HIT
As a convenience for users simply wishing to use some basic functionality of SpecBit, DecayBit and PrecisionBit from the command line, we also provide an additional driver program 3-BIT-HIT. This program doubles as an example of how the three modules can be used in standalone mode, i.e. without the GAMBIT Core or the other modules.
3-BIT-HIT takes a single MSSM model as input via a minimal YAML file 3bithit.in, in either the weak-scale MSSM20atQ or GUT-scale NUHM2 parameterisation. The user can also specify SM parameters, following SLHA2 conventions. The program then solves the RGEs and evaluates pole masses, calculates decay widths and branching fractions of all SM particles and their superpartners, computes electroweak precision observables, and finally outputs the results as a single SLHA file The Higgs-sector decays are computed by DecayBit using FeynHiggs, SM decays are based on native DecayBit functions, and sparticle decays come from SDECAY, via SUSY-HIT and DecayBit. 3-BIT-HIT computes the SUSY contribution to the anomalous magnetic moment of the muon, along with its uncertainty, using calls to GM2Calc from PrecisionBit. All other electroweak precision observables (EWPO) come from FeynHiggs via PrecisionBit, where errors are also assigned to them.
3-BIT-HIT performs a similar function to SUSY-HIT, but computes EWPO in addition to spectra and decays. It also employs FlexibleSUSY instead of SuSpect for spectrum generation, FeynHiggs for Higgs decays rather than HDECAY, and our patched version of SDECAY (see 3.1.3 for details) for sparticle decays. Its source can be found in DecayBit/examples/3bithit.cpp, and it can be built with make 3bithit

Summary
In this paper we have introduced the GAMBIT modules SpecBit, DecayBit and PrecisionBit. These are designed to flexibly integrate publicly available programs for spectrum generation, the calculation of decay widths and branching ratios, and additional precision calculations, like the anomalous magnetic moment of the muon. Together, these modules provide a powerful way to synergise spectrum generators, decay codes and additional precision calculations, allowing users to extract the information produced by them in a common format. We have illustrated this use with the example standalone program 3-BIT-HIT (Sec. 5), which provides a single interface to FlexibleSUSY, FeynHiggs, SDECAY and GM2Calc. At the same time, these modules play a crucial role within the full framework of GAMBIT, where they provide essential information to other packages and important components of the log-likelihood functions used to drive global fits.
Sven Heynemeyer and Sebastian Paßehr for doing so on Feyn-Higgs. PA would like to thank his FlexibleSUSY collaborators Dominik Stöckinger, Jae-hyeon Park, Dylan Harries and especially Alexander Voigt, for many helpful discussions. We warmly thank the Casa Matemáticas Oaxaca, affiliated with the Banff International Research Station, for hospitality whilst part of this work was completed, and the staff at Cyfronet, for their always helpful supercomputing support.

A Physics Models
Here we give a compact physical definition of each of the models currently supported by SpecBit, DecayBit and PrecisionBit, complementing the slightly more codefocussed description in the main GAMBIT paper [10]. GAMBIT is primarily designed for global fits to extensions of the Standard Model (SM) of particle physics. Many calculations in these extensions require SM quantities, and many such extensions make minimal modifications to the SM; examples are the scalar singlet dark matter model (Sec. A.3) and effective weakly-interacting massive particle (WIMP) dark matter (see Ref. [12] for an example in GAMBIT). We therefore begin with the SM itself.

A.1: Standard Model
The SM describes the interactions of all observed fundamental particles. These interactions are invariant under the symmetries of the gauge group, The Lagrangian density is given by, where u Ri , d Ri and e Ri represent the right-handed uptype quark, down-type quark and lepton. Similarly, represent the Higgs, left-handed quark and left-handed lepton doublets respectively. Here i and j specify the generation of the fields, a and b are SU (2) L indices, and ab is the anti-symmetric tensor. The B µν , W A µν and G B µν are field strength tensors for the B µ , W A µ and G B µ fields respectively, which are respectively associated with the U (1) Y , SU (2) L and SU (3) C gauge groups. The kinetic terms for the matter fields include / D = D µ γ µ , where D µ is the usual covariant derivative, which, for example, takes the following form for all fields transforming non-trivially under both SU (2) L and SU (3) C : Here g 1 and Q Y are the GUT-normalised gauge coupling and GUT-normalised charges 14 of U (1) Y . Similarly g 2 and T A are the SU (2) L gauge couplings and generators respectively, and g 3 and λ B are the SU (3) C gauge couplings and generators respectively. Electro-weak symmetry is broken when the Higgs field develops a VEV, v, so that field may be rewritten as where h 0 is the physical Higgs boson and G 0 is the unphysical Goldstone boson that corresponds to the longitudinal mode for the massive Z 0 boson.
14 Hypercharge is often written without GUT normalisation, such that the charge for a field φ is Y φ = 5/3Q φ Y and will be written in the covariant derivative with the gauge coupling g = 3/5g1.

The SM VEV gives masses to the SM particles. The Z and W gauge bosons obtain tree-level masses
The weak mixing angle θ W is given by though it is important to note that the last equality depends on the tree-level relations given in Eq. A.7, so this definition holds only within the renormalisation scheme in which the gauge couplings are defined. An alternative definition can be made in terms of pole masses, Finally, one should note that neither of these correspond exactly to the leptonic effective weak mixing angle, which is an important precision observable. This quantity is relevant to PrecisionBit and is discussed in Sec. 4.2.2.

A.2: Minimal Supersymmetric Standard Model
The MSSM superpotential is given by, where a hat indicates a superfield and i, j specify the generation, while a, b are SU (2) L indices. This includes the Higgsino mixing parameter µ, which appears as a dimension one coefficient of the up-and down-type Higgs superfields, The chiral superfields have dimensionless up-and downtype Yukawa couplings, Y u and Y d with the Higgs superfields and the chiral superfields for the right-handed quarks. Specifically left-handed quarks,Q, interact with both the chiral superfields for the right-handed up quarks,Û c , and the up-type Higgs superfield through Y u and with down quarks,D c , and the down-type Higgs superfields through Y d . Similarly the chiral superfield for left-handed leptonsL has Yukawa interactions, Y e with the chiral superfields for the right-handed leptons, E c andĤ d .
Terms which break supersymmetry without introducing quadratic divergences are given in the soft SUSY breaking Lagrangian, and (m 2 e ) ij respectively. Finally we have soft trilinears for the up-type squarks and Higgs, (T u ) ij , down-type squarks and Higgs, (T d ) ij , and the sleptons and down-type Higgs, (T e ) ij and also the terms that are are soft when there is no singlet in the model, (C u ) ij , (C d ) ij and (C e ) ij . The latter are normally assumed to be zero, even though the MSSM has no scalar singlets as it is difficult to generate these. The soft trilinears may also be rewritten in terms of the following, In the MSSM the neutral components of the up-and down-type scalar Higgs fields, H 0 u and H 0 d respectively, develop VEVs, 14) and the bilinear soft term, b, mixes them, so that they must be rotated to obtain the physical mass eigenstates, where the mixing matrix is commonly expressed in terms of a mixing angle α, Similarly the pseudoscalar Higgs state, A, is formed from the imaginary parts of the Higgs scalar fields where G 0 is the neutral Goldstone boson. At tree level it is straightforward to show that the mass eigenstates can be obtained with a rotation, giving a massless Goldstone boson and an expression for the pseudoscalar mass, The charged Higgs mass eigenstate, H + and charged Goldstone boson, G + are related to the charged scalar Higgs fields through, (A.20) The squarks are rotated into the mass eigenstates as, 21) and the sleptons and sneutrinos are rotated as, The neutralinos are mass eigenstates from combinations of the BinoB 0 , neutral WinoW 3 and the up-and down-type Higgsinos,H 0 u andH 0 d , where we follow the SLHA [1] convention and require the mixing matrix, Z N , to be real, which means the masses of the neutralinos may be negative. The mass matrix of the charginos is diagonalised through a biunitary transformation, such that the two component chargino mass eigenstates are given by,

A.3: Scalar Singlet Dark Matter Model
The scalar singlet model is a minimalistic extension of the SM, with one new scalar charged only under a Z 2 symmetry. This symmetry guarantees that the scalar is a stable dark matter candidate, and restricts the most general permitted renormalisable Lagrangian density on symmetry arguments to be, where µ 2 S is a mass squared term for the singlet, λ hS the Higgs-portal coupling and λ S the quartic self coupling.
The singlet mass has an additional contribution from EWSB through the Higgs-portal coupling, giving a mass, where v is the SM Higgs VEV in Eq. A.5. As the Z 2 symmetry must remain unbroken for stable dark matter we demand that the scalar field, S, does not obtain a VEV and thus no mixing occurs between the Higgs and the new scalar field. This requirement is satisfied for λ S > 0.

B Spectrum generators as backends: SPheno
One of the key features of GAMBIT is its ability to interface with different external codes to perform relevant calculations for a specific scan. These codes, or backends, are linked with GAMBIT in a non-trivial way in order to access their functions direcly via memory, rather than through a shell. Details and descriptions of the backend system and how to create a new backend can be found in the main GAMBIT paper [10]. Consequently, as with other calculations, the generation of a spectrum can be delegated to a backend library which provides the required capability. At the time of the first release, GAMBIT comes along with Flex-ibleSUSY [18] internally implemented as the de facto spectrum generator for MSSM models, but also allows the use of the backended spectrum generator SPheno [19,20]. In this section we describe how SPheno is added as a backend to GAMBIT and how it provides the spectrum capability unimproved_MSSM_spectrum, to serve as a guide for the addition of other spectrum generators in the future, e.g. SOFTSUSY [30][31][32][33] or SuSpect [34]. It is worth noting that SPheno is written in Fortran, so some of the details of the backending process will differ significantly from software written in other languages; SOFTSUSY, for example, is written in C++, so it would need to draw on the classloading tool BOSS (Backend-On-a-Stick Script, described in [10]).

B.1: Installation and import of variables/functions from the backend
The general guideline to create a GAMBIT backend from an external tool, as described in the main GAMBIT paper [10], is to compile it into a shared library in order to provide access to the internal variables and functions at runtime. This thus requires patching the makefile provided by SPheno to create this shared library. Effectively one needs to add the flag -shared to the rule lib/libspheno.so, as well as the flag -fPIC to the compilation rule for every source file. Once compiled, this makes the variables and functions inside SPheno available to use in GAMBIT. However, the main Fortran program SPheno does not get a shared library symbol, and thus all the functions and subroutines inside, such as the one that calculates the spectrum CalculateSpectrum, are not made available to GAMBIT. This is easly fixed by a further patch into the main file SPheno3.f90, to change the Program directive for a Module directive, and commenting all the executed lines in the program, but leaving the subroutines unscathed.
Therefore, all the variables and functions can be accessed by the backend system. These are imported by using the macros BE_VARIABLE and BE_FUNCTION, described in the Backends section of the GAMBIT main paper [10], with the library symbols ''__module_MOD_var'' for a variable var in a module module. Note that Fortran is case insensitive, but all library symbols are written in small caps (except the MOD tag). The files SpecBit/include/gambit/SpecBit/frontends/SPheno.hpp and SpecBit/src/frontends/SPheno.cpp are then created to import these variables and functions, and to give any specifications that are necessary.
Along with the imported variables and functions, a set of convenience functions is created, by the macro BE_CONV_FUNCTION, see Ref. [10], specified in the frontend file SPheno.cpp. These are run_SPheno, ReadingData, InitializeStandardModel, Spectrum_Out and ErrorHandling. In particular the convenience function run_SPheno will effectively be the main function of the backend and will cover the executed lines commented out by the patch, avoiding those that perform input and ouput operations, which will be described in the section below.
Lastly, every backend gets a backend initialisation function, where any overall initialisations for the backend are performed, inside the namespace BE_INI_FUNCTION in the frontend source file. The only initialisation required by SPheno at this stage is the information about the model which is done by the function ModelInUse.

B.2: Input and output, warnings and errors
SPheno follows the SLHA conventions [1,2] for the format of the spectra it generates, and as such all its input and (successful) output happens in the form of files using the SLHA format. Reading and writing from files is inefficient and hence, for the most part, GAMBIT avoids performing any of these tasks at the point level, for it would slow down the scans significantly. Therefore, all the input required and output provided by SPheno is stored and obtained directly to and from the backend variables, as described above, via two convenience functions ReadingData and Spectrum_Out.
The convenience function ReadingData initialises all required SPheno variables for the relevant model and scan. Firstly all the SM parametes are initialised in the convenience function InitializeStandardModel using the struct SMInputs, followed by some internal initialization of variables via the backended functions InitializeLoopFunctions and Set_All_Parameters_0. In addition all the run options, as described in detail in Sec. 2.2.6, are parsed at this step.
Secondly, all the information pertaining to the SLHA blocks MINPAR, EXTPAR, MSL2, MSE2, MSQ2, MSU2, MSD2, TUIN, TDIN and TUIN is pulled from the model parameters and stored in the corresponding SPheno variables. As expected, not every variable is required to be filled for every model, e.g. for the CMSSM, only the variables relating to the MINPAR block are filled, which are obtained from the model parameters M0, M12, TanBeta, SignMu and A0. Variables relating to other SLHA input blocks are either filled somewhere else, SMINPUTS in InitializeStandardModel, and MODSEL in the backend initialization function, as described above; or are ignored for they have no internal use as an input, such as SPINFO.
The output convenience function Spectrum_Out deals with the creation of a Spectrum object from the spectrum generated by SPheno. As was mentioned in Sec. 2.2.4, the backended SPheno version does not provide postgeneration RGE running, like FlexibleSUSY does. Rather, the spectrum is stored in a SLHAstruct variable and then transformed into a static Spectrum object via the function spectrum_from_SLHAea. By default the scale dependent output (SLHA blocks MSOFT, GAUGE, etc) is given at the SUSY scale.
Warnings and error messages in SPheno are also written into a file and, likewise, GAMBIT diverts that output to avoid slowing down the scans. Internally this is done by initializing the variable ErrCan to 0, hereby forcing SPheno to write any output into an unused buffer and not to a file. Error tracking is hence done by use of the integer variable kont, whose value depends on the specific error that occurred. For any non-zero value of kont, the convenience function ErrorHandling appends the corresponding message to the logger and raises an invalid_point exception. Specific details on which values of kont correspond to which error messages can be found in Appendix C of the SPheno manual [20].

B.3: Calculation of the spectrum
The spectrum is calculated in SPheno by using the backended function CalculateSpectrum, provided all required variables are initialised previously. After this function is run, all the output variables are filled (including the error tracking variable kont), and they are written into a Spectrum object, as described above.
Once the spectrum is calculated, it is returned to the function that had this backend requirement, by providing the capability SPheno_MSSMspectrum. The built-in function with this backend requirement is get_MSSM_spectrum_SPheno. This function provide the capability unimproved_MSSM_spectrum, effectively passing along the created Spectrum object to the next level in the dependency tree. The definition of the function get_MSSM_spectrum_SPheno is: where Finputs is a data structure containing the input parameters and the run options.
All of the additional calculations available in the out-of-the-box version of SPheno, such as low energy observables, branching ratios and cross sections, are disabled in this backended version, as their functionality is already covered by other GAMBIT modules: FlavBit [14], DecayBit and ColliderBit [13].

C List of SpecBit SubSpectrumContents definitions
Spectrum information that has been stored in SpecBit may be accessed via the SubSpectrum interface class. The stored information is given in Tables 16 (for low energy SM information), 17 (for the MSSM) and 18 (for the scalar singlet model). To extract this information SpecBit has accessors that have one of the following forms, subspectrum.get(tag,label) subspectrum.get(tag,label,index 1) subspectrum.get(tag,label,index 1,index 2) The arguments label and tag are given in Tables 16, 17 and 18 in columns 3 and 4 respectively. The fifth column in these tables specifies the number of indices that are required as additional arguments, which will be either 0, 1 or 2. When necessary the last two columns in these tables provide the range each index runs over. These indices specify the element in a vector or matrix of parameters that are structured this way.
If the using namespace Par directive is not in use in the scope, then tag must be prefixed with the namespace qualifer Par, e.g. Par::tag.
The SubSpectrum classes provide interfaces to specific backend spectrum generators, while SubSpectrumContents classes define what parameters a SubSpectrum wrapper must provide. In Table 12 we list the SubSpectrumContents definitions which ship with the SpecBit module in GAM-BIT 1.0.0, and describe each of them in detail in the table referenced beside each list item.

D List of SpecBit SubSpectrum wrappers
Here we provide a list of all the SubSpectrum wrapper classes which ship with the GAMBIT 1.0.0 version

SMHiggs SMSimpleSpec
General container for SM spectrum data. Can be constructed from an SLHAea object

SM_slha ScalarSingletDM SimpleSpec
Simple wrapper for scalar singlet DM parameters, plus SM gauge and Yukawa couplings SingletDM   Table 13: List of available "simple" SubSpectrum wrappers, which wrap only parameter data and provide no RGE facilities, and are therefore not connected to any backend spectrum generators. Spectrum data is either entered directly as input parameters or retrieved from an external source, such as an SLHA file. The final column 'Contents' gives the SubSpectrumContents definition (see Appendix C) to which the wrapper conforms.  of SpecBit. Wrappers based on simple parameter containers with no RGE facilities are listed in Table 13, while wrappers that connect to full spectrum calculators are listed in Table 14. The interface to the wrapped content of these classes is defined by their associated SubSpectrumContents classes, which are given in the "Contents" column of each table. The SubSpectrumContents classes are themselves described in Appendix C.

E SubSpectrum wrapper GetterMaps and SetterMaps function signatures
In this appendix we give details of the allowed function pointer signatures for the GetterMaps and SetterMaps objects which consistute the link between SubSpectrum interface functions, and backend spectrum data (see Sec. F.2). In principle the system can be extended to allow any function signature. However, it cannot be done automatically, and thus at present only a limited set of options is available. If the object to which you want to interface does not conform to one of these signatures, you will need to write 'helper' functions in your wrapper which alter the function signature. Some of the function signatures listed are designed to support such helper functions.
The list is given in Table 15. So, for example, to store a function pointer with the signature double f(Model&) in a GetterMaps collection, one needs to put it into the map0 map, as it shown in the example in Sec. F.3. The map0W (see Sec. F.2), map0_extraM and map0_extraI maps work similarly, with function signatures as shown in the table.
Functions associated with indices require an extra step; they must be associated with lists of the allowed indices. This is made easy via the FInfo helper classes, which are listed in table 15. For example, say we wish to attached two-index functions with the signature double Model::f(int,int) to the get interface. The fill_getter_maps function for a wrapper Wrap would then be: where std::set<int> i012 is initialised with a helper function initSet (from the header gambit/ Utils/util_functions.hpp). This set is then used to tell the GetterMaps that each index of the function can accept the values 0, 1 or 2. At present, it is not possible to restrict the allowed index values pair-wise (that is, GAMBIT cannot automatically provide safeguards against users asking for specific invalid combinations of indices).
For more examples of how to add functions to GetterMaps and SetterMaps please see the source code for the various SubSpectrum wrappers that ship with SpecBit. These are listed in Appendix D.

F Adding support for new spectrum calculators
While the GAMBIT modules described in this paper are initially distributed with the MSSM and scalar singlet models already implemented, the great advantage of this framework is the manner in which new models can be systematically added. In this section we describe how to do this. The most simple way to extend SpecBit is to just add a new SubSpectrumContents definition, which is described in Sec. F.1. This creates a new standardised set of string label/tag/index sets to which new SubSpectrum wrappers may choose to conform. This is the first step that is required for adding a wrapper for an entirely new model, but it does not actually provide any real capabilities, it only standardises the interface for wrappers which provide access to the same basic physics model.
The real work of interfacing to new spectrum information is done by writing a SubSpectrum wrapper, which is the standardised GAMBIT interface to any container of spectrum information. This information may be supplied using a new spectrum generator created with FlexibleSUSY, a new model in SARAH / SPheno, or any available public or private spectrum generator for that model (that is written in C/C++ or Fortran; codes written in other languages can at present only be launched via calls to the operating system).
The new spectrum generator then needs to be interfaced with SpecBit. If the spectrum generator is written in C++ and has getters and setters to access the spectrum information (and a run method) then this task can be greatly simplified. It is also possible to use the spectrum objects as a simple container for parameters, with no RGE running facilities (i.e. no actual spectrum generator "hooked up"), which is the simplest use case.
There are therefore a number of different ways that the Spectrum interface class can be connected to parameter data. We will tackle these various cases in order of increasing complexity throughout the following sections.
As well as designing the wrapper class, in order to actually use the new spectrum object one must write GAMBIT module functions that construct the objects and return them as GAMBIT capabilities (running the backend spectrum generator if required). This task will be discussed in Sec. F.7.
As a quick-reference for writers of new wrapper classes, in Sec. 2.4 we provide a checklist of tasks to be completed when writing a wrapper, with references to the subsection of this guide in which the relevant detailed instructions can be found.

F.1: SubSpectrumContents definitions
In this section we describe the mechanism which enforces the consistency of SubSpectrum objects wrapping different spectrum calculators. Recall that SubSpectrum objects are simply a virtual interface class; the real work of interfacing to external spectrum calculators is performed by wrapper classes which derive from SubSpectrum. In the process of defining this wrapper class, one must associate the wrapper with a SubSpectrumContents class. The SubSpectrumContents classes are simple structs that define which parameters are supposed to be contained by wrappers that interface to a calculator for a specific physics model. So, for example, two wrapper classes which wrap two different external spectrum calculators, but which calculate the spectrum for the same physics model, should both be associated with the same SubSpectrumContents class. The shared SubSpectrumContents class then acts as a promise that identical 'get' calls made to each of the two wrappers should retrieve equivalent information. Furthermore, if a wrapper fails to make certain parameters accessible that are defined as part of the SpectrumContents, then a runtime error will occur when the spectrum is constructed. The error will explain that the wrapper does not conform to the declared contents and so cannot be used. This message should serve only to help wrapper writers not to forget any contents; it will not occur when using completed wrappers.
The first release of SpecBit defines the set of SubSpectrumContents classes listed below: and descriptions of each of them can be found in Appendix C. The above code shows the entire declaration for these classes, as can be found in Models/include/gambit/Models/SpectrumContents/Regi-steredSpectra.hpp, with the exception that each class is derived from a base class SubSpectrumContents (omitted here for brevity). The only definition required for these classes is their constructor, which for the above classes are defined in source files with matching names in the Models/src/SpectrumContents/ directory. Defining the constructor, which defines the required parameters  for that spectrum type, is very simple. It is easiest to explain with an example, so below we show the definition for the ScalarSingletDM class: SpectrumContents:: ScalarSingletDM::ScalarSingletDM() { using namespace Par; setName("ScalarSingletDM"); addParameter(mass1, "vev" ); addParameter(dimensionless, "lambda_hS"); addParameter(Pole_Mass, "h0"); addParameter(Pole_Mass, "S" ); } First, a string name for the spectrum contents should be declared via the setName function. After that one simply declares all the parameters which should exist in the SubSpectrum wrapper by specifying the tag/string name/indices by which that parameter should be accessed. In the example the index definition argument is omitted, and so the parameter is assumed to be scalarvalued and require no index. To define an index requirement, one supplies as a final argument a vector of integers, specifying the dimension size of each required index. For example, to define a vector parameter with six entries, and a matrix parameter requiring two indices each of size three, one would write e.g. std::vector<int> v6; v6.push_back(6); std::vector<int> m3x3; m3x3.push_back(3); m3x3.push_back(3); addParameter(Pole_Mass, "∼d" , v6); addParameter(mass2 , "mq2", m3x3); The SubSpectrumContents constructor definitions can always be consulted to check exactly what content is required to be retrievable from any wrapper that conforms to it. For reference purposes we provide tables in Appendix C which describe each of the pre-defined SubSpectrumContents sets, with cross-referencing to theoretical descriptions of each parameter. The descriptions of pre-defined SubSpectrum wrappers in Appendix D also refer to these tables to describe what they contain.

F.2: Wrapping a simple parameter collection
The most basic kind of object for which one may wish to constuct an interface is a simple class that contains parameter values. Such an object is far simpler than the typical case that the SubSpectrum interface is designed to wrap, and so many of its features are not needed. The functionality provided by wrapping these simple objects in a SubSpectrum could be entirely replaced by the standard GAMBIT Model system (see Ref. [10]), which already provides a way to deal with simple parameter containers, however it can be useful to interact with parameters via the Spectrum interface in cases where they might alternatively be provided by a true spectrum generator. It is also useful to examine this case simply to demonstrate the most basic requirements that a wrapper class must fulfill. So let us consider the following simple class, and construct a wrapper which links it to the SubSpectrum interface. This wrapper is in fact implemented in GAM-BIT and provides some scalar-singlet dark matter pa- Weak mixing angle "sinW2" Pole_Mixing 0 n/a n/a αe Fine structure constant "alpha" dimensionless 0 n/a n/a αs Strong coupling constant "alphaS" dimensionless 0 n/a n/a mγ Photon M S mass "gamma" mass1 0 n/a n/a mg Gluon M S mass "g" mass1 0 n/a n/a mu Up-type quark M S masses "u"/"ubar" mass1 1 1,2,3 * n/a m d Down-type quark M S masses "d"/"dbar" mass1 1 1,2,3 * n/a m l Charged lepton M S masses "e-"/"e+" mass1 1 1,2,3 * n/a m h 0 Higgs boson pole mass "h0" Pole_Mass 0 n/a n/a v in (A.5) Higgs VEV "vev" mass1 0 n/a n/a The namespace used here is not important, we specify it just to match the actual code.
To begin wrapping this struct, one must define a specialisation of the SpecTraits template class. The purpose of this class is to communicate essential type information to the SubSpectrum base classes. Suppose that our Bino soft mass "M1" mass1 n/a n/a M2 in (A.12a) Wino soft mass "M2" mass1 n/a n/a M3 in (A.12a) Gluino soft mass "M3" mass1 n/a n/a M u in (A.10) Superpotential Higgs bilinear "Mu" mass1 n/a n/a vu in (A.14) Up-type Higgs VEV "vu" mass1 n/a n/a vd in (A.14) Down-type Higgs VEV "vd" mass1 n/a n/a  Higgs pole mass "h0" Pole_Mass 0 n/a n/a v in (A.5) Higgs VEV "vev" mass1 0 n/a n/a λ in (A.2) Higgs quartic coupling "lambda_h" dimensionless 0 n/a n/a λs in (A.26) Scalar quartic coupling "lambda_S" dimensionless 0 n/a n/a λ hs in (A.26) Higgs quartic coupling "lambda_hS" dimensionless 0 n/a n/a g1 in (A.4) U (1) Y gauge coupling "g1" dimensionless 0 n/a n/a g2 in (A.4) SU (2) W gauge coupling "g2" dimensionless 0 n/a n/a g3 in (A.4) SU (3) C gauge coupling "g3" dimensionless 0 n/a n/a sin 2 θ W in Eq A.8 weak mixing angle "sinW2" dimensionless 0 n/a n/a  wrapper class is to be named ScalarSingletDMSimpleSpec; the associated traits class should then be: Here the namespace is important; the original SpecTraits template is declared in the Gambit namespace, so the specialisation must also live in that namespace. Note that we forward declare the class ScalarSingletDMSimpleSpec, which will be our wrapper, because we need it as the template parameter for the traits class.
The required members of this SpecTraits specialisation are as follows: name() A function returning a std::string name for the wrapper. This is used in error messages so the class name is generally the logical choice. Contents A typedef which identifies the contents definition to which this wrapper will conform (see SubSpectrumContents in sec. F.1) Other SpecTraits members can be defined, however for this simple example we do not need them (the defaults inherited from DefaultTraits will suffice). We will return to this in the more complicated examples of sections F.3 and F.4.
We are now ready to define the wrapper itself. The declaration is as follows: Let us discuss what is going on here. First, the wrapper class must inherit from the Spec class, and provide its own type as the template parameter. This is because we employ the CRTP (curiously-recurring template pattern) for static polymorphism, to allow the base class access to the wrapper member functions.
Second, the wrapper contains an instance of the SingletDMModel object to which we want to interface. This is not necessary, but it is helpful for maintaining encapsulation.
Next is the constructor. The wrapper writer is quite free to do what they like with this; here we simply use it to initialise the member object.
Following the constructor are a series of "getter" and "setter" functions, which retrieve and set the parameter values we are interested in. These are the functions that we will "hook up" to the SubSpectrum "get" and "set" functions. These functions can access the hosted SingletDMModel in our example, and so could also perform extra tasks like calling functions of the SingletDMModel (if it had any) or performing unit conversions. Our later examples will demonstrate tasks like this. For now, it will suffice for these functions to have the following sort of definition: We will refrain from listing the rest of the definitions in this example because they follow a similar pattern. Finally, we get to the key part of the wrapper; the fill_getter_maps and fill_setter_maps functions. These are what define the actual interface to the hosted Model (and Input) objects. They are simply functions that fill a series of map (i.e. std::map) structures in the base class, where these maps are what define which functions are 'hooked up' to the get, set and has functions of the SubSpectrum interface (see Sec. 2.3.3). The operation of these 'filler' functions is most easily seen by looking at the definition of one, so let us examine the fill_getter_maps definition for our example: Here GetterMaps is another type inherited from the base Spec class, and is the main container object that we use to associate a tag/string pair with a function pointer. So the first entry is what makes it possible for a call subspectrum.get(Par::mass1,"vev") to the interface class to in-turn call the get_HiggsVEV member function of Model. The fill_getter_maps() function is used by the base Spec class to initialise a member variable of type GetterMaps, so here we are in fact defining how this member variable will be 'filled'. The fill_setter_maps function works in direct analogy to the fill_getter_maps function, so we will skip discussion of it.
Note that neither of these filler-functions must be defined, the wrapper will compile without them and the interface will simply not accept any tag/string pairs to the get/set functions. For example one may fill the 'getter' maps, but not the 'setter' maps, and it will then simply be impossible to change the parameters in the underlying Model object via the 'set' functions, which is not necessarily problematic. However, if one fails to fill the 'getter' maps in accordance with the declared Contents of the wrapper, then a runtime error will be thrown as soon as the wrapper object constructor is called, which will present a message telling the user that the wrapper is not correctly defined.
There is a subtlety here; note the map0W data member which is accessed from the map_collection variable. This member is what defines the function signature of the function pointer to be stored in the map_collection. This function signature must be known 'in advance' in order for the base class to call the function pointer, so the function pointer needs to be placed in the correct location within the map_collection. In our case the map0W data member expects functions with the signature double Self::function(void); that is, they should be a member function of the wrapper, take no arguments, and return a double. The equivalent map0W member within SetterMaps expects the signature void Self::function(double); that is, functions should be members of the wrapper, should accept a double, and return void.
Functions with a variety of other signatures can also be used, however they must be placed in the correct data member of GetterMaps/SetterMaps in order to work. A full list of the allowed function signatures, the way they should be stored in the GetterMaps/SetterMaps, and the way they can be accessed via the SubSpectrum get and set functions, is given in Appendix E. We will see some of these other signatures in use in the more complicated wrapper examples. When one wants to create a SubSpectrum wrapper for preexisting classes, it may be the case that the pre-existing class already has "getter" and "setter" functions, such that one doesn't need to write new ones in the wrapper. Unfortunately, it is not possible to store functions of arbitrary signature in the GetterMaps and SetterMaps, and the pre-enabled set of permitted signatures is not very large. More can be added, however it is quite technically involved. If you have a special need to add more then please contact the authors for advice. The practical use-cases of directly using external member functions is therefore very limited, however in SpecBit we make use of this feature for our most complicated wrapper; the MSSMSpec wrapper which interfaces to FlexibleSUSY model objects. Furthermore, new autogenerated FlexibleSUSY spectrum generators will come with the required set of getters/setters with the correct function signature, so it is useful to describe this 'shortcut' interface method just for the sake of the FlexibleSUSY case.
Suppose that the external class to be wrapped has the following form: It does not matter where the parameters themselves live at present, we just need the functions that set and retrieve them. In this example we will also show how to deal with functions requiring indices.
These functions can be inserted into the GetterMaps and SetterMaps of the wrapper via the filler functions as follows: where these function signatures match the map0 and map2 members of the GetterMaps and SetterMaps classes, as shown in Appendix E. This appendix also gives an explanation of the FInfo2 helper class and how it is used to help specify the allowed indices to get and set functions.
In order to call these functions via function pointers, the SubSpectrum wrapper needs to know the C++ type of the hosted class. This can be supplied via a typedef in the SpecTraits struct; The wrapper base classes will then be able to use the MyClass type via the Model typedef. In addition, the wrapper needs to have access to an instance of the MyClass type. This must be provided by overloading a special inherited member function get_Model in the wrapper, e.g. Note that here Model is the typedef for MyClass, which is learned via the SpecTraits struct. We assume in this example that an instance of model is carried as a data member of the wrapper class, which is a good idea for maintaining encapsulation, but it is not strictly required. All that is required is that get_Model() return an instance of the Model type. Note also that both const and nonconst versions of get_Model() are required, because the wrapper base classes need both in order to maintain const-correctness when const wrapper objects are used.

F.4: Interfacing with non-class functions
In the previous SubSpectrum wrapper examples, we covered interfacing with member functions of wrapper objects, and with member functions of particular external classes. The former can be designed to access external functions or classes in arbitrary ways and so are the most generally useful, while the latter have limited use due to restrictions on the allowed function signature but are convenient for dealing with certain special cases such as FlexibleSUSY classes. Now, we will deal with interfacing to plain functions that are not members of any class.
As with the other permitted kinds of functions, these are also restricted to specific function signatures. In fact, they do not provide any additional functionality over simply using member functions of the wrapper classes, and wrapper member functions are usually a better choice since they have easier access to data members of the wrapper class. However, they can be useful to avoid code repetition if the target functions are to be used in more than one wrapper, for example if they perform some common unit conversions or simple calculations.
The allowed non-class function signatures are described in Appendix E. They may be connected to GetterMaps and SetterMaps as in the following example: Here, Input is imagined to be a class containing input information used to setup the wrapper, but which one may also wish to access via the SubSpectrum interface along with the rest of the spectrum data. But it is simply an arbitrary class which can be used to pass information to the interface functions.

F.5: Index offsets
In general, SubSpectrum wrapper interfaces use a onebased indexing system. For example when retrieving neutralino pole masses from the MSSM spectrum wrappers the lightest neutralino is retrieved with the index 1. However, not all external functions will follow this convention, so it is useful to have a system for converting between index systems. The wrapper system provides this functionality via the index_offset wrapper class member function. This function should be overridden in the wrapper in order to define an offset to be added to all index values input by users to the get and set functions of that wrapper. Its function signature is: static int index_offset() {return offset;} So, for example, if offset is -1, then when a user calls e.g. subspec.get(Pole_Mass,"∼e",1) the offset will be added to the input index of 1, resulting in 0, before being passed to the external function (which will therefore see the index as 0). If the index_offset function is not overriden then a default null offset of 0 will be applied.
Note that this function should have public access. If the external spectrum generator is a C++ code, then connecting the SubSpectrum wrapper to the RGE running facilities should be fairly easy. There are two main functions that need to be connected: GetScale, and RunToScaleOverride. These functions are not required by the SubSpectrumContents definitions, so whether or not a particular SubSpectrum can perform RGE running will need to be inferred from the GAMBIT capability it is given via whatever module function provides it.
The GetScale function is the simplest, and should return the scale at which all running parameters are defined. If there is an analogous function or parameter defined in the host Model object then the GetScale definition will simply be something like this: Of course the details for hooking up these functions correctly are entirely dependent on how the external code works. The onus is therefore on the wrapper writer to understand the external code. However, the power of this framework is that once the wrapper is defined, writers of other module functions can interact with spectrum information in a totally generic way, even including RGE running. F.7: Constructing and returning Spectrum objects from module functions So far, we have described how to design a new wrapper class for external spectrum data. However, in GAM-BIT, a Spectrum object needs to be constructed from the wrapper and returned via a module function before the spectrum data can be used by other module functions.
Here we discuss the general requirements of this process. The details for constructing the wrapper depend on how the wrapper and its constructor are defined, which are up to the wrapper writer. We will therefore use the SingletDM example wrapper from Sec. F.2 to demonstrate the more general requirements for creating and returning the abstract interface. Let us first write out an example module function which does this, and then we will examine it line-by-line.
First, the module function rollcall declaration: Let us look first at the module function definition, and then afterwards we will examine how the rollcall declaration ensures that the necessary dependencies are made available. The first line gives the function signature, which, as always, returns void and takes one argument by reference, which is the pre-allocated memory for the storage of the function result. Here we will create a Spectrum object so the 'result' type is Spectrum. Next, a dependency on an SMInputs object (as described in Sec. 2.3.5) is retrieved.
Lines 7 to 11 show the construction of an instance of the SingletDMModel class, which is the object that we want to wrap. Here we set the parameters contained in the class by extracting them from the Params map supplied by GAMBIT via the core hierarchical model system (see Sec. 5 of Ref. [10]). On line 14 we then use the SingletDMModel instance to construct our instance of the ScalarSingletDMSimpleSpec wrapper class.
Nearing the end of the module function, we construct on line 16 the full Spectrum object out of the components we have prepared, namely the scalar singlet SubSpectrum object and an SMInputs object. More generally one can also supply the 'low-energy' SubSpectrum component, however this shortened constructor will automatically generate an SMSimpleSpec (see Appendix D) wrapper based on SMINPUTS to fulfil this role. The final argument to the Spectrum constructor is the set of model parameters which were available to this Spectrum when it was constructed, which are retrieved via the Params map.
Finally, the newly-constructed object is moved into the result space of the module function, from where it can be distributed to other module functions by the GAMBIT dependency resolver. Now, let us return to the module function rollcall declaration and see how it requests the dependencies that we require. The first two lines simply declare the name for the capability that this function provides. The next two (define/start function) declare the C++ name of the function (which of course must match the function as defined elsewhere). Note here that we also declare the result type for this function as Spectrum.
Next, we declare the dependency we require. Here it is declared as SMINPUTS, with type SMInputs. This will then be made available to our module function by the GAMBIT core. The next few lines declare a joint dependence on model parameters from both StandardModel_Higgs and SingletDM, that is, they declare that this module function requires parameters from both of these models. For more details on this syntax please see Ref. [10].
With the rollcall declaration done we are finished, and the Spectrum object created by the module function get_SingDM_spec is made available to the rest of GAMBIT under the capability name SingletDM_spectrum.

F.8: Controlling wrapper lifetimes
That covers the construction of a simple wrapper, however there are a few more subtleties regarding the lifetime (in the sense of time between construction/destruction of the C++ objects as the code runs) of the various wrapper objects to be discussed. There are several constructors for the Spectrum object, and it is important to choose the correct one depending on how you want the member SubSpectrum objects to be treated. The options are listed below, with the argument accepting the GAMBIT parameter container replaced with <Param> for brevity. Spectrum(SubSpectrum const le, SubSpectrum const he, const SMInputs& smi, const std) -Construct new object, wrapping existing SubSpectrum objects. If the original objects are prematurely destructed then attempting to access them via the Spectrum interface will cause a segmentation fault.

G Worked example of writing a SubSpectrum wrapper
In this section we provide a stripped-down tutorial-style example of how to define a new SubSpectrum wrapper and return a Spectrum object from a module function.
Comments and instructions will be kept to a bare minimum, and instead we will simply refer to the appropriate sections of the paper body in which further details can be found.

G.1: FlexibleSUSY MSSM wrapper
The first step in adding a new FlexibleSUSY spectrum generator to SpecBit is to run FlexibleSUSY and generate the C++ code for the new spectrum generator. It then needs to be added to the GAM-BIT cmake build system, and the appropriate headers included in certain key SpecBit headers. As mentioned in Sec. 2.1.2 we provide instructions for these steps as separate documentation, which can be found in gambit/doc/Adding_FlexibleSUSY_Models.txt, and they will change when future technical improvements to the GAMBIT BOSS [10] system permit FlexibleSUSY to be automatically backended. Thus we begin this example by assuming that these steps have been completed. From this point we can begin following the checklist given in Sec 2.4. 1) Choose or create a SubSpectrumContents definition (Sec. F.1) -The purpose of the SubSpectrum classes is to enforce consistency between wrappers which wrap the same essential physics, so one should always try to use an existing definition if possible. However, for the purposes of this example we will pretend that the definition for the MSSM did not previously exist. Thus, in gambit/Models/include/gambit/Models/ SpectrumContents/RegisteredSpectra.hpp we add a line of code: = initVector(4); std::vector<int> v6 = initVector(6); std::vector<int> m2x2 = initVector(2,2); std::vector<int> m3x3 = initVector(3,3); std::vector<int> m4x4 = initVector(4,4); std::vector<int> m6x6 = initVector(6,6); // Spectrum parameters // tag, name, shape addParameter(mass2, "BMu" , scalar); addParameter(mass2, "mHd2", scalar); addParameter(mass2, "mHu2", scalar); The SpectrumContents that we will use are now defined.
2) Write a SpecTraits specialization -In this example, the getter and setter functions are defined in an external object, so the example in Sec. F.4 is relevant. It is best to define our SpecTraits class in the header file for the SubSpectrum wrapper class, because the definition needs to be available to the C++ compiler when when actual SubSpectrum class is defined. We thus begin a new header  We assume here that the new FlexibleSUSY spectrum generator has been defined such that the soft masses are input at a user-specifed scale Q, such that the matching GAMBIT model parameters are MSSM63atQ. We therefore declare ALLOW_MODELS(MSSM63atQ), so that we will be able to access these parameters in our module function. We will also need SM SLHA2 parameters, but it is more convenient to obtain them via an SMInputs object rather than directly from GAMBIT model parameters, so we declare a dependency on SMINPUTS, which can be provided by an existing SpecBit module function.

H Glossary
Here we explain some terms that have specific technical definitions in GAMBIT.
backend An external code containing useful functions (or variables) that one might wish to call (or read-/write) from a module function . backend function A function contained in a backend. It calculates a specific quantity indicated by its capability . Its capability and call signature are defined in the backend's frontend header . backend requirement A declaration that a given module function needs to be able to call a backend function or use a backend variable , identified according to its capability and type(s). Backend requirements are declared in module functions' entries in rollcall headers. backend variable A global variable contained in a backend . It corresponds to a specific quantity indicated by its capability. Its capability and type are defined in the backend's frontend header . BOSS The Backend-On-a-Stick script, used for preprocessing C++ backend code to allow GAMBIT to dynamically load classes from it. capability A name describing the actual quantity that is calculated by a module or backend function. This is one possible place for units to be noted; the other is in the documented description of the capability (see Sec. 10.7 of Ref. [10]). dependency A declaration that a given module function needs to be able to access the result of another module function, identified according to its capability and type. Dependencies are declared in module functions' entries in rollcall headers. dependency resolver The component of the GAMBIT Core that performs dependency resolution . dependency resolution The process by which GAM-BIT determines the module functions, backend functions and backend variables needed and allowed for a given scan, connects them to each others' dependencies and backend requirements, and determines the order in which they must be called. dependency tree A result of dependency resolution ; a directed acyclic graph of module functions connected by resolved dependencies. See Fig. 5 of Ref. [10] for an example.
frontend The interface between GAMBIT and a given backend , consisting of a frontend header plus optional source files and type headers. frontend header The C++ header in which the frontend to a given backend is declared. module A subset of GAMBIT functions following a common theme, able to be compiled into a standalone library. Although module often gets used as shorthand for physics module , this term technically also includes the GAMBIT scanning module ScannerBit. module function A function contained in a physics module . It calculates a specific quantity indicated by its capability and type , as declared in the module's rollcall header . It takes only one argument, by reference (the quantity to be calculated), and has a void return type. physics module Any module other than ScannerBit, containing a collection of module functions following a common physics theme. rollcall header The C++ header in which a given physics module and its module functions are declared. type A general fundamental or derived C++ type, often referring to the type of the capability of a module function .