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

TL;DR

This paper introduces GAMBIT modules SpecBit, DecayBit, and PrecisionBit, which enable seamless integration of spectrum generators, decay codes, and precision observable calculations for comprehensive particle physics analyses.

Contribution

The paper presents a modular framework that links various computational tools for particle physics, allowing flexible, consistent, and automated calculations within GAMBIT.

Findings

Modules enable automatic combination of existing codes

Results are presented in a unified, accessible format

The framework supports comprehensive global fits in particle physics

Abstract

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.

Figures

Tables

Table 1. Table 1 : SM cap abilities provided by SpecBit . The cap ability SM_spectrum is understood to provide running parameters in the M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} scheme, along with pole masses.

Capability	Function (Return Type): Brief Description	Dependencies
SMINPUTS	get_SMINPUTS (SMInputs): Provide Standard Model parameters in SLHA2 input conventions.
qedqcd_subspectrum	get_QedQcd_spectrum (const SubSpectrum*): Create QedQcdWrapper version of SubSpectrum from SMInputs structures.	SMINPUTS
SM_spectrum	get_SM_spectrum (Spectrum): Create Spectrum object from QedQcdWrapper and SM Higgs parameters.	SMINPUTS
Higgs_Couplings	SM_higgs_couplings (HiggsCouplingsTable): Construct a table of SM Higgs couplings.	Higgs_decay_rates

Table 2. Table 2 : MSSM cap abilities provided by SpecBit . The cap abilities MSSM_spectrum and unimproved_MSSM_spectrum are understood to provide MSSM couplings and soft masses in the D ​ R ¯ ¯ 𝐷 𝑅 \overline{DR} scheme, along with pole masses, with low-energy Standard Model inputs defined in the M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} scheme. Note also that module functions exist in PrecisionBit with cap ability MSSM_spectrum , which depend on the unimproved_MSSM_spectrum ; see Table 10 and Sec.  4.2.3 for more details.

Capability	Function (Return Type): Brief Description	Dependencies	Backend requirements
unimproved_MSSM	get_CMSSM_spectrum(Spectrum): Make an MSSM spectrum object with FlexibleSUSY using CMSSM boundary conditions.	SMINPUTS	FlexibleSUSY
_spectrum
	get_MSSMatMGUT_spectrum(Spectrum): Make an MSSM spectrum object with FlexibleSUSY using soft breaking masses input at the GUT scale.	SMINPUTS	FlexibleSUSY
	get_MSSMatQ_spectrum (Spectrum): Make an MSSM spectrum object with FlexibleSUSY from soft breaking masses input at (user-defined) scale $Q$.	SMINPUTS	FlexibleSUSY
	get_MSSM_spectrum_SPheno(Spectrum): Make an MSSM spectrum object with SPheno for any supported MSSM model.	SMINPUTS	SPheno
	get_MSSM_spectrum_from_SLHAfile (Spectrum): Construct an MSSM spectrum from an SLHA file. Wraps it up in MSSMSimpleSpec; i.e. no RGE running possible. Designed for testing against benchmark points, but may be a useful last resort for interacting with “difficult” spectrum generators.
	get_unimproved_MSSM_spectrum_as_SLHAea (SLHAstruct): Convert an unimproved MSSM spectrum into an SLHAstruct object.	unimproved_MSSM_spectrum
	get_unimproved_MSSM_spectrum_as_map (map_str_dbl): Convert an unimproved MSSM spectrum into a std::map, for saving via the GAMBIT printer system.	unimproved_MSSM_spectrum
MSSM_spectrum	get_MSSM_spectrum_as_map (map_str_dbl): Convert an MSSM spectrum into a std::map, for saving via the GAMBIT printer system.	MSSM_spectrum
	get_MSSM_spectrum_as_SLHAea (SLHAstruct): Turn an MSSM spectrum into an SLHAstruct object.	MSSM_spectrum
SM_subspectrum	get_SM_SubSpectrum_from_MSSM_Spectrum (const SubSpectrum*): Extract the “low-energy” SubSpectrum from an MSSM Spectrum object.	MSSM_spectrum
SMlike_Higgs	most_SMlike_Higgs_MSSM (int): Return the PDG code of the most SM-like of the neutral Higgs bosons in the MSSM.	MSSM_spectrum
_PDG_code
Higgs_Couplings	MSSM_higgs_couplings_pwid (HiggsCouplingsTable): Construct a table of Higgs couplings from Higgs partial decay widths.	MSSM_spectrum
		SMlike_Higgs_PDG_code
		Reference_SM_Higgs_decay_rates
		Reference_SM_other_Higgs_decay
		_rates
		Reference_SM_A0_decay_rates
		Higgs_decay_rates
		h0_2_decay_rates
		A0_decay_rates
		H_plus_decay_rates
		t_decay_rates
	MSSM_higgs_couplings_FH (HiggsCouplingsTable): Construct a table of Higgs couplings from couplings calculated by FeynHiggs.	MSSM_spectrum
		SMlike_Higgs_PDG_code
		Reference_SM_Higgs_decay_rates
		Reference_SM_other_Higgs_decay
		_rates
		Reference_SM_A0_decay_rates
		Higgs_decay_rates
		h0_2_decay_rates
		A0_decay_rates
		H_plus_decay_rates
		t_decay_rates
		FH_Couplings_output
FH_Couplings	FH_Couplings (fh_Couplings): Higgs couplings, as computed by FeynHiggs	prec_HiggsMasses	FeynHiggs
_output
prec_HiggsMasses	FH_HiggsMasses (fh_HiggsMassObs): Higgs masses and mixings with theoretical uncertainties, as computed by FeynHiggs.	FH_MSSMMasses	FeynHiggs
FH_MSSMMasses	FH_MSSMMasses (fh_MSSMMassObs): SUSY masses and mixings, as computed by FeynHiggs.		FeynHiggs

Table 3. Table 3 : Capabilities provided by SpecBit for the SingletDM model. The cap ability SingletDM_spectrum is understood to provide running parameters in the M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} scheme, along with pole masses, with low-energy Standard Model inputs also defined in the M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} scheme.

Capability	Function (Return Type): Brief Description	Dependencies	Options
			(Type)
SingletDM_spectrum	get_SingletDM_spectrum_simple (Spectrum): Create simple container Spectrum object from SMInputs structures, SM Higgs parameters, and the SingletDM parameters.	SMINPUTS
	get_SingletDM_spectrum_FS (Spectrum): Create Spectrum object with RGE facilities for the scalar singlet model using FlexibleSUSY.	SMINPUTS	FlexibleSUSY
	get_SingletDM_spectrum_as_map (map_str_dbl): Convert a SingletDM spectrum into a std::map for printing.	SingletDM_spectrum
Higgs_Couplings	SingletDM_higgs_couplings_pwid (HiggsCouplingsTable): Construct a table of Higgs couplings from Higgs partial decay widths.	SingletDM_spectrum
		Reference_SM_Higgs
		_decay_rates
		Higgs_decay_rates
vacuum_stability	find_min_lambda (triplet<double>): Run the quartic Higgs coupling to high scales, returning a triplet containing the expected lifetime of the electroweak vacuum, the instability scale, and the perturbativity of the couplings up to the requested scale.	SingletDM_spectrum	set_high_scale
			(double)
			check_perturb_pts
			(int)
check_perturb_min	get_check_perturb_min_lambda (double): Check that dimensionless couplings are sufficiently small up to the smaller of the Planck scale and the scale at which the quartic Higgs coupling is minimised and non-positive.	vacuum_stability
_lambda
VS_likelihood	get_likelihood (double): Get the a priori likelihood of no electroweak vacuum decays having occurred in the past light cone.	vacuum_stability
expected_lifetime	get_expected_lifetime (double): Return the expected lifetime of the electroweak vacuum in units of years.	vacuum_stability

Table 4. Table 4 : SM decay functions provided by DecayBit . None of these functions has any backend requirements. Here X is either empty or equal to bar , Y ∈ { \in\{ plus , minus } } \} and Z ∈ { \in\{ 0 , plus , minus } } \} . All functions with X = bar , Y = minus or Z = minus have a dependency on their C ​ P 𝐶 𝑃 CP -conjugate equivalent, with X empty, Y = plus or Z = plus .

Capability	Function (Return Type): Brief Description	Dependencies	Options (Type)
tX_decay_rates	tX_decays(DecayTable::Entry): Computes all SM decays of the $t$/$\overline{t}$ quark.
mu_Y_decay_rates	mu_Y_decays(DecayTable::Entry): Computes all SM decays of the ${\mu }^{+}$/${\mu }^{-}$ lepton.
tau_Y_decay_rates	tau_Y_decays(DecayTable::Entry): Computes all SM decays of the ${\tau }^{+}$/${\tau }^{-}$ lepton.
pi_Z_decay_rates	pi_Z_decays(DecayTable::Entry): Computes all SM decays of the ${\pi }^0$, ${\pi }^{+}$ and ${\pi }^{-}$ mesons.
rho_Z_decay_rates	rho_Z_decays(DecayTable::Entry): Computes all SM decays of the ${\rho }^0$, ${\rho }^{+}$ and ${\rho }^{-}$ mesons.
eta_decay_rates	eta_decays(DecayTable::Entry): Computes all SM decays of the $\eta$ meson.
omega_decay_rates	ompega_decays(DecayTable::Entry): Computes all SM decays of the $\omega$ meson.
Z_decay_rates	Z_decays(DecayTable::Entry): Computes all SM decays of the $Z$ boson.
W_Y_decay_rates	W_Y_decays(DecayTable::Entry): Computes all SM decays of the $W^{+}$/$W^{-}$ gauge boson.
Reference_SM	Ref_SM_Higgs_decays_table (DecayTable::Entry): Computes all decays of the Higgs boson in the SM, by interpolating in the tables of Ref. YellowBook13 .	mh	higgs_minmass(double)
_Higgs_decay_rates			higgs_maxmass(double)
	Ref_SM_Higgs_decays_FH (DecayTable::Entry): Computes all decays of the Higgs boson in the SM, using Higgs couplings obtained from FeynHiggs.	SMlike_Higgs_PDG_code
		FH_Couplings_output
		SLHA_pseudonyms
Higgs_decay_rates	SM_Higgs_decays(DecayTable::Entry): Returns all decays of the SM Higgs boson.	Reference_SM
		_Higgs_decay_rates

Table 5. Table 5 : Scalar singlet functions provided by DecayBit . Neither of these functions has any backend requirements nor options.

Capability	Function (Return Type): Brief Description	Dependencies
Higgs_decay_rates	SingletDM_Higgs_decays(DecayTable::Entry): Computes all decays of the Higgs boson, including to scalar singlets.	SingletDM_spectrum
		Reference_SM_Higgs_decay_rates
lnL_Higgs_invWidth	lnL_Higgs_invWidth_SMlike(double): Computes the log-likelihood of the Higgs invisible width.	Higgs_decay_rates

Table 6. Table 6 : MSSM decay functions provided by DecayBit . None have any options of note. Here A and B are mass indices that run from 1 to 2 and 1 to 4 respectively, C is a helicity label ( l or r ), and X and Y are as in Table  4 (i.e. empty/ bar and plus / minus ).

Capability	Function (Return Type): Brief Description	Dependencies	Backend requirements
t_decay_rates	FH_t_decays(DecayTable::Entry): Computes MSSM decays of the $t$ quark.	FH_Couplings_output
Higgs_decay_rates	MSSM_h0_1_decays(DecayTable::Entry): MSSM decays of the lightest $CP$-even Higgs with SUSY-HIT.	SLHA_pseudonyms	SUSY-HIT
	FH_MSSM_h0_1_decays(DecayTable::Entry): MSSM decays of the lightest $CP$-even Higgs with FeynHiggs.	SLHA_pseudonyms
		FH_Couplings_output
h0_2_decay_rates	h0_2_decays(DecayTable::Entry): Decays of the heaviest $CP$-even MSSM Higgs with SUSY-HIT.	SLHA_pseudonyms	SUSY-HIT
	FH_h0_2_decays(DecayTable::Entry): Decays of the heaviest $CP$-even MSSM Higgs with FeynHiggs.	SLHA_pseudonyms
		FH_Couplings_output
A0_decay_rates	A0_decays(DecayTable::Entry): Decays of the $CP$-odd MSSM Higgs boson with SUSY-HIT.	SLHA_pseudonyms	SUSY-HIT
	FH_A0_decays(DecayTable::Entry): Decays of the $CP$-odd MSSM Higgs boson with FeynHiggs.	SLHA_pseudonyms
		FH_Couplings_output
HY_decay_rates	HY_decays: Decays of the MSSM charged Higgs with SUSY-HIT.	SLHA_pseudonyms	SUSY-HIT
	FH_HY_decays: Decays of the MSSM charged Higgs with FeynHiggs.	SLHA_pseudonyms
		FH_Couplings_output
Reference_SM_other	Ref_SM_other_Higgs_decays_table (DecayTable::Entry): Computes all decays of an SM Higgs with mass equal to the mass of the least SM-like $CP$-even Higgs boson in the MSSM, by interpolating in the tables of Ref. YellowBook13 and extensions computed with HDECAY 6.51.	SMlike_Higgs_PDG_code
_Higgs_decay_rates		MSSM_spectrum
	Ref_SM_other_Higgs_decays_FH(DecayTable::Entry): Computes all decays of an SM Higgs with mass equal to the mass of the least SM-like $CP$-even Higgs boson in the MSSM, using Higgs couplings obtained from FeynHiggs.	SMlike_Higgs_PDG_code
		FH_Couplings_output
		SLHA_pseudonyms
Reference_SM_A0	Ref_SM_A0_decays_table(DecayTable::Entry): Computes all decays of an SM Higgs with mass equal to the mass of the $CP$-odd Higgs boson in the MSSM, by interpolating in the tables of Ref. YellowBook13 and extensions computed with HDECAY 6.51.	SMlike_Higgs_PDG_code
_decay_rates		MSSM_spectrum
	Ref_SM_A0_decays_FH(DecayTable::Entry): Computes all decays of an SM Higgs with mass equal to the mass of the $CP$-odd Higgs boson in the MSSM, using Higgs couplings obtained from FeynHiggs.	SMlike_Higgs_PDG_code
		FH_Couplings_output
		SLHA_pseudonyms
gluino_decay_rates	gluino_decays(DecayTable::Entry): Computes decays of the gluino.	SLHA_pseudonyms	SUSY-HIT
charginoY_A_decay_rates	charginoY_A_decays(DecayTable::Entry): Computes decays of the charginos.	SLHA_pseudonyms	SUSY-HIT
neutralino_B_decay_rates	neutralino_B_decays(DecayTable::Entry): Computes decays of the neutralinos.	SLHA_pseudonyms	SUSY-HIT
supX_C_decay_rates	supX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{u}$ squarks.		SUSY-HIT
sdownX_C_decay_rates	sdownX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{d}$ squarks.		SUSY-HIT
scharmX_C_decay_rates	scharmX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{c}$ squarks.		SUSY-HIT
sstrangeX_C_decay_rates	sstrangeX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{s}$ squarks.		SUSY-HIT
stopX_A_decay_rates	stopX_A_decays(DecayTable::Entry): Computes decays of the $\tilde{t}$ squarks.	SLHA_pseudonyms	SUSY-HIT
sbottomX_A_decay_rates	sbottomX_A_decays(DecayTable::Entry): Computes decays of the $\tilde{b}$ squarks.	SLHA_pseudonyms	SUSY-HIT
selectronX_C_decay_rates	selectronX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{e}$ sleptons.		SUSY-HIT
smuonX_C_decay_rates	smuonX_C_decays(DecayTable::Entry): Computes decays of the $\tilde{\mu }$ sleptons.		SUSY-HIT
stauX_A_decay_rates	stauX_A_decays(DecayTable::Entry): Computes decays of the $\tilde{\tau }$ sleptons.	SLHA_pseudonyms	SUSY-HIT
snuX_electronl_decay_rates	snuX_electronl_decays(DecayTable::Entry): Computes decays of the ${\tilde{\nu }}_{e_L}$ sneutrinos.		SUSY-HIT
snuX_muonl_decay_rates	snuX_muonl_decays(DecayTable::Entry): Computes decays of the ${\tilde{\nu }}_{{\mu }_L}$ sneutrinos.		SUSY-HIT
snuX_taul_decay_rates	snuX_taul_decays(DecayTable::Entry): Computes decays of the ${\tilde{\nu }}_{{\tau }_L}$ sneutrinos.	SLHA_pseudonyms	SUSY-HIT

Table 7. Table 7 : Collector and helper module functions provided by DecayBit . Here X refers to all particles that exist in a given model. Note that str is a GAMBIT -specific typedef of std :: string .

Capability	Function (Return Type): Brief Description	Dependencies	Options (Type)
decay_rates	all_decays(DecayTable): Computes all decays of all particles in a given model.	X_decay_rates	drop_SLHA_file(bool)
			SLHA_output_filename(str)
	all_decays_from_SLHA(DecayTable): Reads decays from SLHA files.		SLHA_decay_filenames
			(std::vector<str>)
SLHA1_violation	check_first_sec_gen_mixing(int): Checks if first and second generation sfermions mix at more than a given %.	MSSM_spectrum	gauge_mixing_tolerance(double)
SLHA_pseudonyms	get_mass_es_pseudonyms (DecayBit::mass_es_pseudonyms): Computes nearest-match SLHA2 mass eigenstates for each SLHA1 sfermion weak/family eigenstate.	MSSM_spectrum	gauge_mixing_tolerance(double)
			gauge_mixing_tolerance
			_invalidates_point_only(bool)

Table 8. 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 𝑊 W and h ℎ h masses, which depend directly on the calculated W 𝑊 W and h ℎ 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 ColliderBit are generally to be preferred to lnL_h_mass_chi2 , as they simultaneously take account of correlated constraints on the mass and couplings.

Capability	Function (Return Type): Brief Description	Dependencies	Options (Type)
lnL_alpha_em	lnL_alpha_em_chi2(double): Computes the log-likelihood of the fine-structure constant, ${\alpha }_{\mathrm{em}}\left(m_Z\right)$, in the $\overline{MS}$ scheme.	SMINPUTS
lnL_alpha_s	lnL_alpha_s_chi2(double): Computes the log-likelihood of the strong coupling constant, ${\alpha }_s\left(m_Z\right)$, in the $\overline{MS}$ 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 $\overline{MS}$ masses of the $u$, $d$ and $s$ quarks at $\mu =2$ GeV.	SMINPUTS	mud_obs(double)
			mud_obserr(double)
			msud_obs(double)
			msud_obserr(double)
			ms_obs(double)
			ms_obserr(double)
lnL_mcmc	lnL_mcmc_chi2(double): Computes the log-likelihood of the $\overline{MS}$ mass $m_c\left(m_c\right)$.	SMINPUTS
lnL_mbmb	lnL_mbmb_chi2(double): Computes the log-likelihood of the $\overline{MS}$ mass $m_b\left(m_b\right)$.	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 9. Table 9 : Summary of PrecisionBit module functions providing the W 𝑊 W and h ℎ h masses from different GAMBIT   Spectrum inputs. These functions have no options nor backend requirements.

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.	MSSM_spectrum
		SMlike_Higgs_PDG_code

Table 10. 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 .

Capability	Function (Return Type): Brief Description	Dependencies	Backend requirements
FH_Precision	FH_PrecisionObs(fh_PrecisionObs): Calculate precision observables in the MSSM with FeynHiggs.	FH_Couplings_output	FeynHiggs
deltarho	FH_precision_deltarho(double): Retrieve parameter $\Delta \rho$ from FeynHiggs precision calculations.	FH_Precision
prec_mw	FH_precision_mw(double): Retrieve $W$ mass in the MSSM from FeynHiggs precision calculations.	FH_Precision
prec_sinW2_eff	FH_precision_sinW2(double): Retrieve effective leptonic weak mixing angle ${\sin }^2{\theta }_{W,\mathrm{eff}}$ from FeynHiggs precision calculations.	FH_Precision
edm_e	FH_precision_edm_e(double): Retrieve electron electric dipole moment from FeynHiggs precision calculations.	FH_Precision
edm_n	FH_precision_edm_n(double): Retrieve neutron electric dipole moment from FeynHiggs precision calculations.	FH_Precision
edm_hg	FH_precision_edm_hg(double): Retrieve mercury electric dipole moment from FeynHiggs precision calculations.	FH_Precision
MSSM_spectrum	make_MSSM_precision_spectrum_H_W(Spectrum): Function to provide an updated MSSM spectrum with precision $W$ and Higgs masses.	unimproved_MSSM_spectrum
		prec_mw
		prec_HiggsMasses
	make_MSSM_precision_spectrum_W(Spectrum): Function to provide an updated MSSM spectrum with precision $W$ mass.	unimproved_MSSM_spectrum
		prec_mw
	make_MSSM_precision_spectrum_none(Spectrum): Function to present an unimproved MSSM spectrum as if it is a precision spectrum.	unimproved_MSSM_spectrum
muon_gm2	FH_precision_gm2(double): Retrieve MSSM contribution to $g-2$ from FeynHiggs precision calculations.	FH_Precision
	GM2C_SUSY(double): Calculate MSSM contribution to $g-2$ using GM2Calc.	MSSM_spectrum	GM2Calc
	SI_muon_gm2(double): Calculate MSSM contribution to $g-2$ using SuperIso.	SuperIso_modelinfo	SuperIso
SP_PrecisionObs	SP_PrecisionObs(double): Calculate precision observables in the MSSM with SUSY-POPE (experimental).		SUSY-POPE

Table 11. Table 11 : PrecisionBit module functions providing log-likelihood calculations for electroweak precision observables.

Capability	Function (Return Type): Brief Description	Dependencies
lnL_sinW2_eff	lnL_sinW2_eff_chi2(double): Computes the log-likelihood associated with ${\sin }^2{\theta }_{W.\mathrm{eff}}$.	prec_sinW2_eff
lnL_gm2	lnL_gm2_chi2(double): Computes the log-likelihood associated with muon $g-2$.	muon_gm2
		muon_gm2_SM
lnL_deltarho	lnL_deltarho_chi2(double): Computes the log-likelihood associated with $\Delta \rho$.	deltarho
muon_gm2_SM	gm2_SM_ee(triplet<double>): 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 ${\tau }^{+}{\tau }^{-}$ data.

Table 12. Table 12 : List of available SubSpectrumContents definitions.

Name	Plain English Description	Table
SM	SM pole masses and QED$\times$QCD running parameters	16
SM_slha	SM pole masses and SLHA SMINPUTS MSbar masses (light quarks)	16
SMHiggs	SM Higgs sector; pole mass and VEV	16
MSSM	63 parameter general MSSM in SLHA2 conventions	17
SingletDM	Scalar singlet dark matter plus SM Higgs sector	18

Table 13. 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.

Name	Plain English Description	Contents
MSSMSimpleSpec	General container for MSSM spectrum. Can be constructed from an SLHAea object	MSSM
SLHASimpleSpec	Base class for wrappers based on SLHAea.	n/a
SMHiggs  SimpleSpec	Very simple wrapper for SM Higgs parameters	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 14. Table 14 : List of available “full” SubSpectrum wrappers, which connect directly to backend spectrum generator code and provide RGE running facilities. The final column ’Contents’ gives the SubSpectrumContents definition (see Appendix C ) to which the wrapper conforms.

Name	Plain English Description	Contents
MSSMSpec	General templated wrapper for all FlexibleSUSY MSSM models	MSSM
QedQcdWrapper	Wraps QED$\times$QCD calculator from SOFTSUSY	SM
SingletDMSpec	Wraps FlexibleSUSY spectrum for scalar singlet model	SingletDM

Table 15. Table 15 : List of data members of the GetterMaps and SetterMaps classes (which are std :: map s), and the corresponding signatures of function pointers which can be stored in those maps. Note, however, that functions requiring index input cannot be stored directly in the map; they must first be associated with sets of allowed values for those indices using the listed “Helper class”; see text. In the table, the name Self refers to the SubSpectrum wrapper class itself (that is, these maps accept function pointers to members functions of the wrapper class).

Data member name	Type	Helper class	“Getter” function signature	“Setter” function signature
map0	fmap0	none	double Model::f()	void Model::f(double)
map1	fmap1	FInfo1	double Model::f(int)	void Model::f(int,double)
map2	fmap2	FInfo2	double Model::f(int,int)	void Model::f(int,int,double)
map0W	fmap0W	none	double Self::f()	void Self::f(double)
map1W	fmap1W	FInfo1W	double Self::f(int)	void Self::f(double,int)
map2W	fmap2W	FInfo2W	double Self::f(int,int)	void Self::f(double,int,int)
map0_extraM	fmap0_extraM	none	double f(Model&)	void f(Model&,double)
map1_extraM	fmap1_extraM	FInfo1M	double f(Model&,int)	void f(Model&,double,int)
map2_extraM	fmap2_extraM	FInfo2M	double f(Model&,int,int)	void f(Model&,double,int,int)
map0_extraI	fmap0_extraI	none	double f(Input&)	void f(Input&,double)
map1_extraI	fmap1_extraI	FInfo1I	double f(Input&,int)	void f(Input&,double,int)
map2_extraI	fmap2_extraI	FInfo2I	double f(Input&,int,int)	void f(Input&,double,int,int)

Table 16. Table 16 : The spectrum contents available via the get function for the wrappers conforming to the SM and SM_slha (Standard Model) contents definitions (Sec. F.1 ), except the Higgs pole mass and VEV, which are available from wrappers conforming to the SMHiggs definition and which contain only those two parameters. We found this split to be convenient since BSM models commonly alter the Higgs sector. The above quantities may all be obtained via subspectrum . get ( Par :: tag , label , index 1 , index 2 ) where label and tag are given in columns 3 and 4 respectively and subspectrum is defined as in Sec.  2.3.3 . The final two arguments are optional and are filled when necessary by the indices from the allowed range given in columns 6 and 7. ∗ These indices only apply to SM contents; for SM_slha wrappers only u_1 , d_1 , and d_2 are available. Note that the weak mixing angle, provided in the form sin 2 ⁡ θ W pole superscript 2 subscript superscript 𝜃 pole 𝑊 \sin^{2}\theta^{\textrm{pole}}_{W} , corresponds to the pole definition from Eq.  23 . The scale at which all M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} parameters are defined is returned by subspectrum . GetScale () (see Sec.  2.3.3 and Fig.  6 ), unless the wrapper object does not support RGE running, in which case the scales should be in accordance with SLHA2 SMINPUTS definitions (this is the convention followed in the SMSimpleSpec wrapper class, for example. See Appendix  D for full lists of the currently available wrapper classes). Note that if definitions in accordance with SLHA2 SMINPUTS are desired, one can always access the parameters in the SMInputs structure of the parent Spectrum object rather than performing the required RGE running; for example spectrum . get_SMInputs (). mS retrieves the strange quark M ​ S ¯ ¯ 𝑀 𝑆 \overline{MS} mass at 2 GeV. See Fig.  7 for the definition of the SMInputs class, and see Fig.  3 and 4 for the structure of Spectrum objects.

Quantity	Plain English Description	string	tag	\pbox3cmnumber
of indices	\pbox3cmindex 1
range	\pbox3cmindex 2
range
$m_{\gamma }^{pole}$	Photon pole mass	"gamma"	Pole_Mass	0	n/a	n/a
$m_g^{pole}$	Gluon pole mass	"g"	Pole_Mass	0	n/a	n/a
$m_Z^{pole}$	Z pole mass	"Z0"	Pole_Mass	0	n/a	n/a
$m_W^{pole}$	W pole mass	"W+"/"W-"	Pole_Mass	0	n/a	n/a
$m_t^{pole}$	Top quark pole mass	"u_3"/"ubar_3"/ "t"/"tbar"	Pole_Mass	0	n/a	n/a
$m_b^{pole}$	Bottom quark pole mass	"d_3"/"dbar_3"/ "b"/"bbar"	Pole_Mass	0	n/a	n/a
$m_e^{pole}$	Electron pole mass	"e-"/"e+"	Pole_Mass	0	n/a	n/a
$m_{\mu }^{pole}$	Muon pole mass	"mu-"/"mu+"	Pole_Mass	0	n/a	n/a
$m_{\tau }^{pole}$	Tau pole mass	"tau-"/"tau+"	Pole_Mass	0	n/a	n/a
$m_{l_i}^{pole}$	Charged lepton pole masses	"e-"/"e+"	Pole_Mass	1	1,2,3	n/a
$m_{{\nu }_i}^{pole}$	Neutrino pole masses	"nu"/"nubar"	Pole_Mass	1	1,2,3	n/a
${\sin }^2{\theta }_W^{\text{pole}}$ in (23)	Weak mixing angle	"sinW2"	Pole_Mixing	0	n/a	n/a
${\alpha }_e$	Fine structure constant	"alpha"	dimensionless	0	n/a	n/a
${\alpha }_s$	Strong coupling constant	"alphaS"	dimensionless	0	n/a	n/a
$m_{\gamma }$	Photon $\overline{MS}$ mass	"gamma"	mass1	0	n/a	n/a
$m_g$	Gluon $\overline{MS}$ mass	"g"	mass1	0	n/a	n/a
$m_u$	Up-type quark $\overline{MS}$ masses	"u"/"ubar"	mass1	1	1,2,3∗	n/a
$m_d$	Down-type quark $\overline{MS}$ masses	"d"/"dbar"	mass1	1	1,2,3∗	n/a
$m_l$	Charged lepton $\overline{MS}$ 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 (19)	Higgs VEV	"vev"	mass1	0	n/a	n/a

Table 17. Table 17 : The spectrum contents available via the get function for wrappers conforming to the MSSM spectrum contents definition (Sec. F.1 ). These quantities may be obtained via subspectrum . get ( Par :: tag , label , index 1 , index 2 ) where label and tag are given in columns 3 and 4 respectively and subspectrum is defined as in Sec.  2.3.3 . The final two arguments are optional and are filled when necessary by the indices from the allowed range given in the last two columns. Note that weak mixing angle here is defined in the D ​ R ¯ ¯ 𝐷 𝑅 \overline{DR} scheme and is simply reconstructed from the gauge couplings using the expression given in Eq.  22 .

Quantity	Plain English Description	string	tag	\pbox3cmindex 1
range	\pbox3cmindex 2
range
$m_{\tilde{g}}^{pole}$	Gluino pole masses	"~g"	Pole_Mass	n/a	n/a
$m_{W^{\pm }}^{pole}$	W boson pole mass	"W+"	Pole_Mass	n/a	n/a
$m_{A^0}^{pole}$	Psedoscalar Higgs pole mass	"A0"	Pole_Mass	n/a	n/a
$m_{H^{+}}^{pole}$	Charged Higgs pole masses	"H+"	Pole_Mass	n/a	n/a
$m_{h^0}^{pole}$	$CP$-even Higgs pole masses	"h0"	Pole_Mass	{1,2}	n/a
$m_{\tilde{u}}^{pole}$	up-type squark pole masses	"~u"	Pole_Mass	{1,2,3,4,5,6}	n/a
$m_{\tilde{d}}^{pole}$	down-type squark pole masses	"~d"	Pole_Mass	{1,2,3,4,5,6}	n/a
$m_{\tilde{e}}^{pole}$	slepton pole masses	"~e"	Pole_Mass	{1,2,3,4,5,6}	n/a
$m_{\tilde{\nu }}^{pole}$	sneutrino pole mass	"~nu"	Pole_Mass	{1,2,3}	n/a
$m_{{\chi }^{\pm }}^{pole}$	chargino pole masses	"~chi+"	Pole_Mass	{1,2}	n/a
$m_{{\chi }^0}^{pole}$	neutralino pole masses	"~chi0"	Pole_Mass	{1,2,3,4}	n/a
$Z_{\tilde{u}}$ in (35)	up-type squark mixing	"~u"	Pole_Mixing	{1,2,3,4,5,6}	{1,2,3,4,5,6}
$Z_{\tilde{d}}$ in (35)	down-type squark mixing	"~d"	Pole_Mixing	{1,2,3,4,5,6}	{1,2,3,4,5,6}
$Z_{\tilde{e}}$ in (36)	slepton mixing	"~e"	Pole_Mixing	{1,2,3,4,5,6}	{1,2,3,4,5,6}
$Z_{\tilde{\nu }}$ in (36)	sneutrino mixing	"~nu"	Pole_Mixing	{1,2,3}	{1,2,3}
$Z_N$ in (37)	neutralino mixing	"~chi0"	Pole_Mixing	{1,2,3,4}	{1,2,3,4}
$U_{+}$ in (38)	chargino mixing	"~chi+"	Pole_Mixing	{1,2}	{1,2}
$U_{-}$ in (39)	chargino mixing	"~chi-"	Pole_Mixing	{1,2}	{1,2}
$Z_H$ in (29)	CP even Higgs mixing	"h0"	Pole_Mixing	{1,2}	{1,2}
$Z_A$ in 31	CP odd Higgs mixing	"A0"	Pole_Mixing	{1,2}	{1,2}
$B\mu$ in (26c)	soft Higgs bilinear	"BMu"	mass2	n/a	n/a
$m_{H_u}^2$ in (26c)	soft Higgs mass squared	"mHu2"	mass2	n/a	n/a
$m_{H_d}^2$ in (26c)	soft Higgs mass squared	"mHd2"	mass2	n/a	n/a
$m_A^2$ in (33)	$\overline{DR}$ pseudoscalar mass squared	"mA2"	mass2	n/a	n/a
${𝐦}_Q^{\mathrm{𝟐}}$ in (26d)	soft squark mass squared	"mq2"	mass2	{1,2,3}	{1,2,3}
${𝐦}_u^{\mathrm{𝟐}}$ in (26d)	soft up-squark mass squared	"mu2"	mass2	{1,2,3}	{1,2,3}
${𝐦}_d^{\mathrm{𝟐}}$ in (26d)	soft down-squark mass squared	"md2"	mass2	{1,2,3}	{1,2,3}
${𝐦}_L^{\mathrm{𝟐}}$ in (26d)	soft slepton mass squared	"ml2"	mass2	{1,2,3}	{1,2,3}
${𝐦}_e^{\mathrm{𝟐}}$ in (26d)	soft slepton mass squared	"me2"	mass2	{1,2,3}	{1,2,3}
$M_1$ in (26a)	Bino soft mass	"M1"	mass1	n/a	n/a
$M_2$ in (26a)	Wino soft mass	"M2"	mass1	n/a	n/a
$M_3$ in (26a)	Gluino soft mass	"M3"	mass1	n/a	n/a
$Mu$ in (24)	Superpotential Higgs bilinear	"Mu"	mass1	n/a	n/a
$vu$ in (28)	Up-type Higgs VEV	"vu"	mass1	n/a	n/a
$vd$ in (28)	Down-type Higgs VEV	"vd"	mass1	n/a	n/a
${𝐓}_u$ in (26e)	Trilinear	"TYu"	mass1	{1,2,3}	{1,2,3}
${𝐓}_d$ in (26e)	Trilinear	"TYd"	mass1	{1,2,3}	{1,2,3}
${𝐓}_e$ in (26e)	Trilinear	"TYe"	mass1	{1,2,3}	{1,2,3}
$g_1$ in (18)	$U{\left(1\right)}_Y$ gauge coupling222222GUT normalised	"g1"	dimensionless	n/a	n/a
$g_2$ in (18)	$SU{\left(2\right)}_L$ gauge coupling	"g2"	dimensionless	n/a	n/a
$g_3$ in (18)	$SU{\left(3\right)}_C$ gauge coupling	"g3"	dimensionless	n/a	n/a
$\tan \beta$	ratio of VEVs $v_u/v_d$	"tanbeta"	dimensionless	n/a	n/a
${\sin }^2{\theta }_W$ in (22)	weak mixing angle	"sinW2"	dimensionless	n/a	n/a
${𝐘}_u$ in (24)	Up-type quark Yukawa	"Yu"	dimensionless	{1,2,3}	{1,2,3}
${𝐘}_d$ in (24)	Down-type quark Yukawa	"Yd"	dimensionless	{1,2,3}	{1,2,3}
${𝐘}_e$ in (24)	lepton Yukawa	"Ye"	dimensionless	{1,2,3}	{1,2,3}

Table 18. Table 18 : The spectrum contents available via the get function for wrappers conforming to the SingletDM (scalar singlet) spectrum contents definition. These quantities may be obtained via subspectrum . get ( Par :: tag , label , index 1 , index 2 ) where label and tag are given in columns 3 and 4 respectively and subspectrum is defined as in Sec.  2.3.3 . The final two arguments are optional and are filled when necessary by the indices from the allowed range given in columns 6 and 7.

Quantity	Plain English Description	string	tag	\pbox3cmnumber
of indices	\pbox3cmindex 1
range	\pbox3cmindex 2
range
$m_S^{pole}$	Scalar pole mass	"s0"	Pole_Mass	0	n/a	n/a
$m_{h^0}^{pole}$	Higgs pole mass	"h0"	Pole_Mass	0	n/a	n/a
$v$ in (19)	Higgs VEV	"vev"	mass1	0	n/a	n/a
$\lambda$ in (16)	Higgs quartic coupling	"lambda_h"	dimensionless	0	n/a	n/a
${\lambda }_s$ in (40)	Scalar quartic coupling	"lambda_S"	dimensionless	0	n/a	n/a
${\lambda }_{hs}$ in (40)	Higgs quartic coupling	"lambda_hS"	dimensionless	0	n/a	n/a
$g_1$ in (18)	$U{\left(1\right)}_Y$ gauge coupling	"g1"	dimensionless	0	n/a	n/a
$g_2$ in (18)	$SU{\left(2\right)}_W$ gauge coupling	"g2"	dimensionless	0	n/a	n/a
$g_3$ in (18)	$SU{\left(3\right)}_C$ gauge coupling	"g3"	dimensionless	0	n/a	n/a
${\sin }^2{\theta }_W$ in Eq. 22	weak mixing angle	"sinW2"	dimensionless	0	n/a	n/a
${𝐘}_u$ in (16)	Up-type quark Yukawa	"Yu"	dimensionless	2	{1,2,3}	{1,2,3}
${𝐘}_d$ in (16)	Down-type quark Yukawa	"Yd"	dimensionless	2	{1,2,3}	{1,2,3}
${𝐘}_e$ in (16)	lepton Yukawa	"Ye"	dimensionless	2	{1,2,3}	{1,2,3}

Equations

$$V(h)=\frac{\lambda (\mu )}{4}{h}^4$$

$$\Gamma \approx {\Gamma }_0{e}^{-B/\hslash },$$

$$B=S_E\equiv 2{\pi }^2{\int }_0^{\infty }{\rho }^3\mathrm{d}\rho \left[\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}\rho }\right)}^2+V\right]$$

$$\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{\rho }^2}+\frac{3}{\rho }\frac{\mathrm{d}\varphi }{\mathrm{d}\rho }=\frac{\mathrm{d}V}{\mathrm{d}\varphi }.$$

$$h(\rho )=\sqrt{\frac{2}{\left|\lambda \right|}}\frac{2R}{{\rho }^2+R^2}.$$

$$S_E=\frac{8{\pi }^2}{3\left|\lambda \right|},$$

$$\Gamma \approx (T_u{\Lambda }_B{)}^4{e}^{-S_E},$$

$$p\approx {\left(e^{140}\frac{{\Lambda }_B}{M_{\text{Pl}}}\right)}^4\exp \left(-\frac{8{\pi }^2}{3\left|\lambda ({\Lambda }_B)\right|}\right).$$

$$\mathcal{L}=\exp \left[-{\left(e^{140}\frac{{\Lambda }_B}{M_{\text{Pl}}}\right)}^4\exp \left(-\frac{8{\pi }^2}{3\left|\lambda ({\Lambda }_B)\right|}\right)\right].$$

$${\Gamma }_{h\to SS}=\frac{{\lambda }_{hS}^2{v}_0^2}{32\pi m_h}{\left(1-4m_S^2/m_h^2\right)}^{1/2}.$$

$${\Gamma }_{\mathrm{tot}}={\Gamma }_{\mathrm{tree}}(1+{\Delta }_1),$$

$${\Gamma }_{\mathrm{tot}}=\frac{{\Gamma }_{\mathrm{tree}}}{1-{\Delta }_1},$$

$${\mathcal{B}}_X={\Gamma }_{X,\mathrm{tree}}(1+{\Delta }_{X,1})/{\Gamma }_{\mathrm{tot}}.$$

$${\sigma }_{\Delta \rho }^{\mathrm{theo}}=\max ({\sigma }_{m_W}^{\mathrm{theo}}/m_W,{\sigma }_{{\sin }^2{\theta }_{W,\mathrm{eff}}}^{\mathrm{theo}}/{\sin }^2{\theta }_{W,\mathrm{eff}}),$$

$$G_{SM}=SU(3{)}_{\mathrm{C}}\times SU(2{)}_{\mathrm{L}}\times U(1{)}_{\mathrm{Y}}.$$

$$\begin{array}{rl}{\mathcal{L}}_{SM}=& {\sum }_{i=1,3}i({\overline{Q}}_i/D{Q}_i+{\overline{L}}_i/D{L}_i+{\overline{u}}_{Ri}/D{u}_{Ri}\\ & +{\overline{d}}_{Ri}/D{d}_{Ri}+{\overline{e}}_{Ri}/D{e}_{Ri})\\ & +{\sum }_{i,j=1,3}[{\sum }_{a,b=1,2}({\mathbf{Y}}_u{)}_{ij}{ϵ}_{ab}{Q}_i^a({\varphi }^{†}{)}^b{u}_{Rj}\\ & +({\mathbf{Y}}_d{)}_{ij}{\overline{Q}}_i\varphi d_{Rj}+({\mathbf{Y}}_e{)}_{ij}{\overline{L}}_i\varphi e_{Rj}+h.c.]\\ & +|D_{\mu }\varphi {|}^2-{\mu }^2|\varphi {|}^2-\lambda |\varphi {|}^4\\ & -\frac{1}{4}{B}_{\mu \nu }{B}^{\mu \nu }-\frac{1}{4}{W}_{\mu \nu }^A{W}^{\mu \nu A}-\frac{1}{4}{G}_{\mu \nu }^B{G}^{\mu \nu B}\end{array}$$

$$\varphi =\left(\begin{array}{c}{\varphi }^{+}\\ {\varphi }^0\end{array}\right),Q_i=\left(\begin{array}{c}{u}_{Li}\\ d_{Li}\end{array}\right),L_i=\left(\begin{array}{c}{\nu }_{Li}\\ e_{Li}\end{array}\right),$$

$$D_{\mu }={\partial }_{\mu }-i{g}_1{Q}_Y{B}_{\mu }-i{g}_2{W}_{\mu }^A{T}^A-i{g}_3{\lambda }^B{G}_{\mu }^B.$$

$$⟨\varphi ⟩=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ v\end{array}\right),$$

$${\varphi }^0=\frac{1}{\sqrt{2}}(v+h^0+i{G}^0).$$

$$m_Z=\frac{1}{2}\sqrt{(g^{\prime 2}+g_2^2)}v{M}_W=\frac{1}{2}{g}_{2}v.$$

$${\sin }^2{\theta }_W\equiv \frac{g^{\prime 2}}{g^{\prime 2}+g_2^2}=1-\frac{M_W^2}{m_Z^2},$$

$${\sin }^2{\theta }_W^{\text{on-shell}}\equiv 1-\frac{M_W^{\text{(pole)}2}}{m_Z^{\text{(pole)}2}}.$$

$$\begin{gathered} \hat{W}={ϵ}_{ab}\{{\sum }_{i,j=1,3}[({\mathbf{Y}}_u{)}_{ij}{\hat{Q}}_i^a{\hat{H}}_u^b{\hat{U}}_j^c-({\mathbf{Y}}_d{)}_{ij}{\hat{Q}}_i^a{\hat{H}}_d^b{\hat{D}}_j^c \\ -({\mathbf{Y}}_e{)}_{ij}{\hat{L}}_i^a{\hat{H}}_d^b{\hat{E}}_j^c]-\mu {\hat{H}}_u^a{\hat{H}}_d^b\}, \end{gathered}$$

$${\hat{H}}_u=\left(\begin{array}{c}{\hat{H}}_u^{+}\\ {\hat{H}}_u^0\end{array}\right){\hat{H}}_d=\left(\begin{array}{c}{\hat{H}}_d^0\\ {\hat{H}}_d^{-}\end{array}\right).$$

$${\mathcal{L}}_{\mathrm{soft}}=$$

$$-\frac{i}{2}\left[M_1^{\prime }{\overline{\stackrel{~}{B}}}^0{\gamma }_5{\tilde{B}}^0+M_2^{\prime }{\overline{\stackrel{~}{W}}}_A{\gamma }_5{\tilde{W}}_A+M_3^{\prime }{\overline{\stackrel{~}{g}}}_B{\gamma }_5{\tilde{g}}_B\right]$$

$$\displaystyle-\Big{[}b{H}^{a}_{u}{H}_{da}+\mathrm{h.c.}\Big{]}-m^{2}_{H_{u}}|H_{u}|^{2}-m^{2}_{H_{d}}|H_{d}|^{2}$$

$$+{\sum }_{i,j=1,3}\{-[{\tilde{Q}}_i^{†}({\mathbf{m}}_Q^2{)}_{ij}{\tilde{Q}}_j+{\tilde{d}}_{\mathrm{R}i}^{†}({\mathbf{m}}_d^2{)}_{ij}{\tilde{d}}_{\mathrm{R}j}$$

$$\displaystyle+\tilde{u}_{\mathrm{R}i}^{\dagger}(\mathbf{m}^{\mathbf{2}}_{u})_{ij}\tilde{u}_{\mathrm{R}j}+\tilde{L}_{i}^{\dagger}(\mathbf{m}^{\mathbf{2}}_{L})_{ij}\tilde{L}_{j}+\tilde{e}_{\mathrm{R}i}^{\dagger}(\mathbf{m}^{\mathbf{2}}_{e})_{ij}\tilde{e}_{\mathrm{R}j}\Big{]}$$

$$+-{ϵ}_{ab}[({\mathbf{T}}_u{)}_{ij}{\tilde{Q}}_i^a{H}_u^b{\tilde{u}}_{\mathrm{R}j}^{†}-({\mathbf{T}}_d{)}_{ij}{\tilde{Q}}_i^a{H}_d^b{\tilde{d}}_{\mathrm{R}j}^{†}$$

$$-({\mathbf{T}}_e{)}_{ij}{\tilde{L}}_i^a{H}_d^b{\tilde{e}}_{\mathrm{R}j}^{†}+h.c.]$$

$$-{ϵ}_{ab}[({\mathbf{C}}_u{)}_{ij}{\tilde{Q}}_i^a{H}_d^{\ast b}{\tilde{u}}_{\mathrm{R}j}^{†}-({\mathbf{C}}_d{)}_{ij}{\tilde{Q}}_i^a{H}_u^{\ast b}{\tilde{d}}_{\mathrm{R}j}^{†}$$

$$-({\mathbf{C}}_e{)}_{ij}{\tilde{L}}_i^a{H}_u^{\ast b}{\tilde{e}}_{\mathrm{R}j}^{†}+h.c.]\}.$$

Full text

††journal: Eur. Phys. J. C

Contents

@starttoc

toc

1   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 Skands:2003cj and SLHA2 Allanach:2008qq 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 Djouadi:2006bz for the MSSM and NMSSMTools Ellwanger:2004xm ; Ellwanger:2005dv ; Ellwanger:2006rn 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. Allanach:2003jw ; Allanach:2004rh ; Staub:2015aea ; Drechsel:2016htw ).

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 gambit , 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 ScannerBit 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 GAMBIT modules calculate dark matter observables (DarkBit DarkBit ), high-energy collider signatures (ColliderBit ColliderBit ) and quantities from flavour physics (FlavBit FlavBit ). An extended description of the structure, features and abilities of GAMBIT can be found in Ref. gambit .

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 CMSSM ; MSSM ; SSDM have been performed.

The three modules have been initially set up to exploit several external codes111When any of these external codes are used as part of GAMBIT or the SpecBit, DecayBit and PrecisionBit modules, the references for that code listed here should be cited along with this manual. A full list of references associated with all external codes currently utilised as GAMBIT backends can be found in the file README.md in the main directory of the GAMBIT source code.:

• •

FlexibleSUSY Athron:2014yba and SPheno Porod:2003um ; Porod:2011nf for spectrum generation, performing calculations in the dimensional reduction ($\overline{DR}$) scheme.

• •

FeynHiggs Heinemeyer:1998yj ; Heinemeyer:1998np ; Degrassi:2002fi ; Frank:2006yh ; Hahn:2013ria ; Bahl:2016brp for additional Higgs and $W$ mass calculations in a mixed $\overline{DR}/$on-shell (OS) scheme.

• •

HDECAY 111We interface with SDECAY and HDECAY via SUSY-HIT Djouadi:2006bz .Djouadi:1997yw ; Spira:1997dg ; Butterworth:2010ym and FeynHiggs Heinemeyer:1998yj ; Heinemeyer:1998np ; Degrassi:2002fi ; Frank:2006yh ; Hahn:2013ria for Higgs decays.

• •

SDECAY 22footnotemark: 2Muhlleitner:2003vg for sparticle decays

• •

FeynHiggs, SuperIso Mahmoudi:2007vz ; Mahmoudi:2008tp ; Mahmoudi:2009zz and GM2Calc gm2calc for calculation of additional precision observables.

The models and backends supported by SpecBit, DecayBit and PrecisionBit will be continually updated, and are expected to grow rapidly. There are immediate plans to add interfaces to the spectrum generators SOFTSUSY Allanach:2001kg ; Allanach:2009bv ; Allanach:2011de ; Allanach:2014nba and SuSpect Djouadi:2002ze for the MSSM, next-to-minimal SOFTSUSY Allanach:2013kza , NMSSMCALC Ender:2011qh ; Graf:2012hh ; Baglio:2013iia ; King:2015oxa and NMSPEC Ellwanger:2006rn for the NMSSM, and many additional models via FlexibleSUSY Athron:2014yba and SARAH / SPheno Staub:2008uz ; Staub:2010jh ; Staub:2012pb ; Staub:2013tta ; Goodsell:2014bna ; Porod:2011nf . For decays in the MSSM, HFOLD Frisch:2010gw and SFOLD Hlucha:2011yk could also be added. We are also adding SUSYHD arXiv:1504.05200 as a first example for a forthcoming Mathematica 333http://www.wolfram.com/mathematica interface and we will also add similar pure and hybrid effective field theory calculations of HSSUSY Bagnaschi:2015pwa ; Athron:2016fuq and FlexibleEFTHiggs Athron:2016fuq from FlexibleSUSY. NMSSM decays will be obtainable from NMSSMCALC Ender:2011qh ; Graf:2012hh ; Baglio:2013iia ; King:2015oxa and NMHDECAY Ellwanger:2004xm ; Ellwanger:2005dv . Decays in other models can be added from SARAH / SPheno and the upcoming FlexibleDecay extension of FlexibleSUSY.

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 PrecisionBit, 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 license444http://opensource.org/licenses/BSD-3-Clause. Note that fjcore Cacciari:2011ma and some outputs of FlexibleSUSY Athron:2014yba (incorporating routines from SOFTSUSY Allanach:2001kg ) are also shipped with GAMBIT 1.0. These code snippets are distributed under the GNU General Public License (GPL; http://opensource.org/licenses/GPL-3.0), with the special exception, granted to GAMBIT by the authors, that they do not require the rest of GAMBIT to inherit the GPL., and can be downloaded from gambit.hepforge.org.

2   SpecBit

All information in GAMBIT about the spectrum of particle masses and their couplings comes from module functions in SpecBit. This includes the pole masses and mixings of all physical states in the model, scheme dependent parameters, such as those defined in $\overline{DR}$ Siegel:1979wq ; Capper:1979ns ; Jack:1994rk or $\overline{MS}$ Bardeen:1978yd schemes555Note that for a given model, if a user wants to support backend codes that require input parameters in a different scheme to the one implemented in SpecBit, then the user must add a function that converts between the two schemes. Consistency can then be ensured by specifying a scheme-specific spectrum as a dependency. Currently there are no conversions between different schemes in SpecBit, because all backends and native calculations can accept $\overline{DR}$ and $\overline{MS}$ parameters as inputs., Higgs couplings and basic SM inputs (such as the top quark pole mass). 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 $\overline{DR}$ or $\overline{MS}$ 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 666The SM data extractable from SpecBit are shown in Table 16. 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 Allanach:2001kg .

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 $\overline{DR}$ parameters of the MSSM, and an additional scale at which the soft SUSY-breaking $\overline{DR}$ parameters are defined. Note that when there are constraints relating the parameters to each other, varying this scale is not equivalent to a reparametrisation, 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.

SpecBit currently implements several spectrum generators using FlexibleSUSY 1.5.1 (with SARAH 4.9.1). 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. These models represent the most general formulation of the $CP$-conserving MSSM with the couplings given in Eq. 26f set to zero (see Sec. A.2 for the full MSSM Lagrangian). MSSM63atMGUT is defined at the scale where the gauge couplings unify, and MSSM63atQ is defined at a user-specified scale $Q$. These are currently the most general SUSY models of the GAMBIT model hierarchy (which is described in Sec. 5 of Ref. gambit ).

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 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 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 $\overline{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 gambit , and in Appendix B for the specific case of SPheno.

For each of these spectrum generators the EWSB conditions are used to fix $|\mu|$ 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\beta=\frac{v_{u}}{v_{d}}$).

The pole masses and mixings calculated by the spectrum generator are stored internally in the GAMBIT Spectrum wrapper in a format that follows SLHA2 conventions Allanach:2008qq allowing full family mixing for sfermions.

Spectrum generators for the scalar singlet dark matter model

The scalar singlet mass spectrum is available in two forms. One is a via a simple container object that is set up without any spectrum generation, in other words the relevant pole masses and coupling, which govern the new physics effects at the TeV scale, are input directly from scanner-generated parameters. The other option is a spectrum object that interfaces with FlexibleSUSY to calculate the spectrum with full RGE running capabilities.

The simple container spectrum is the most efficient option for a range of calculations in GAMBIT that re- quire only the masses and coupling at a fixed scale. In this case, the only new model parameters required for other module functions are $m_{S}$ and $\lambda_{hs}$. In the container object these parameters are stored directly from the GAMBIT model parameters.

If radiative corrections and/or renormalisation group running of model parameters is necessary, then a fully-calculable 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 $\overline{MS}$ quantities, while all other SM masses are taken according to SLHA2 conventions Allanach:2008qq . 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 $\lambda$ (Eq. 16), 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. gambit for more details), or left as the top mass by default. Optimally, this scale should be set to minimise logarithmic contributions to the RGE $\beta$ functions with the largest coefficient.

** User interface and options**

General settings

In this section we describe how to use SpecBit to run available spectrum generators via a GAMBIT input YAML file.777YAML is the language that GAMBIT input files are written in; see www.yaml.org. At the most rudimentary level, spectrum information can simply be written to disk for one model point in an SLHA2-like format, using SLHAea.888SLHAea is a C++ class allowing internal representation of SLHA files; see fthomas.github.io/slhaea. At a more advanced level, the information can be written to disk for every point in a GAMBIT scan, via the GAMBIT printer system (for details, see Sec. 9 of the GAMBIT Core paper and manual gambit ), and analysed 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 Sec. 2.3. For details on YAML usage of the capabilities and associated options described here, please see the GAMBIT manual gambit 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 gambit . 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 PDB . Users who prefer the $W$ mass returned to respect the tree-level relationship with $m_{Z}$ and $\sin^{2}\theta_{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.

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.

The capability Higgs_Couplings is described in Sec. 2.6.

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.

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 spectra 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 Belanger2013a ; Alanne2014 ; Khan2014 . 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.

• •

check_high_scale: takes a double to define the scale that the couplings are run to (given in units of GeV) after spectrum generation. Default $M_{\text{Pl}}=1.22\times 10^{19}$ GeV.

• •

invalid_point_fatal: kill the whole scan if a FlexibleSUSY error is encountered during spectrum generation. Default false.

• •

FlexibleSUSY options: see Sec. 2.2.5.

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:

• get_MSSMatMGUT_spectrum,
• get_MSSMatQ_spectrum, and
• get_CMSSM_spectrum.

These module functions provide the unimproved_MSSM_spectrum capability in the form of a Spectrum object by calling one of the FlexibleSUSY spectrum generators described in Sec. 2.1.2. The module functions all depend on SMINPUTS provided in the form of an SMInputs structure. These spectra are of the most sophisticated type, providing full two-loop RGE running, as described in Sec. 2.3.3. Options can be supplied to FlexibleSUSY via the input YAML file, and are discussed in Sec. 2.2.5.

Additionally the module function

• get_MSSM_spectrum_SPheno

also provides the capability unimproved_MSSM_spectrum, but in this case obtains the information from a different spectrum generator: SPheno. As with the module functions that call FlexibleSUSY, this module function depends on an SMInputs object and provides a Spectrum object, selecting the specific model and running mode of the backended SPheno at runtime. In contrast to the FlexibleSUSY case, however, this Spectrum object does not provide RGE running. Is is simply a static spectrum equivalent to the contents of an SLHAea object (SLHAstruct), which is a simple container of SLHA-like information (described in more detail in Sec. 2.3.3 under SLHAea output). The options available for the backend version of SPheno used by this module function are in Sec. 2.2.6.

The module function

• get_MSSM_spectrum_from_SLHAfile

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 gambit , 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(\beta-\alpha)>\cos(\beta-\alpha)$, where $\alpha$ is the mixing angle defined in Eq. 30; 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,

• FH_MSSMMasses,
• FH_HiggsMasses and
• FH_Couplings,

respectively provide MSSM masses with capability FH_MSSMMasses, Higgs masses with capability prec_HiggsMasses, and Higgs couplings with capability FH_Couplings_output.

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 Athron:2014yba .

• •

precision_goal: relative error for the calculation of the $\overline{DR}$ spectrum. Default $1\times 10^{-4}$.

• •

max_iterations: maximum number of iterations for the two-scale algorithm. Default $-10\log_{10}$(precision_goal).

• •

calculate_sm_masses: compute SM pole masses during spectrum calculation. Default false.

• •

pole_mass_loop_order: number of loops in self energies for pole mass calculation. Default 2.

• •

ewsb_loop_order: number of loops in radiative corrections to the EWSB conditions in Eq. 33 of Ref. Athron:2014yba . Default 2.

• •

beta_loop_order: loop order for RGEs. Default 2.

• •

threshold_corrections_loop_order: loop order for threshold corrections as defined in Eqs. 6–13 of Ref. Athron:2014yba . Default 2.

• •

Two-loop Higgs pole mass corrections are controlled for specific terms with a boolean input, these are by default all set to true. The available terms are:

• –

use_higgs_2loop_at_as: $\mathcal{O}(\alpha_{t}\alpha_{s})$,

• –

use_higgs_2loop_ab_as: $\mathcal{O}(\alpha_{b}\alpha_{s})$,

• –

use_higgs_2loop_at_at: $\mathcal{O}(\alpha_{t}^{2}+\alpha_{t}\alpha_{b}+\alpha_{b}^{2})$,

• –

use_higgs_2loop_atau_atau: $\mathcal{O}(\alpha_{\tau}^{2})$,

where for $x=t,b,\tau$, $\alpha_{x}=y_{x}^{2}/(4\pi)$ and $y_{x}$ is the corresponding Yukawa coupling, and $\alpha_{s}$ is the strong coupling constant.

Note that the included spectrum generators created by FlexibleSUSY perform a fixed-order calculation of the Higgs mass and we do not currently support additional FlexibleSUSY codes999Such options will be added at a later date, however in the intervening period, an advanced user may also add these codes themselves, following the instructions given in this manual. HSSUSY and FlexibleEFTHiggs that resum large logarithms. However we refer the reader to Ref. Athron:2016fuq for a discussion of how the large logarithms fortuitously cancel in the fixed order FlexibleSUSY calculation, such that the uncertainty is much lower at higher scales in the MSSM than one would naively expect.

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. gambit 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 Porod:2011nf , 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.

• •

SPA_convention: use the SPA convention AguilarSaavedra:2005pw , effectively setting the parameter output scale to 1 TeV. Default false.

• •

GUTScale: fixed value for the GUT scale. Default 0.0, which causes the GUT scale to be obtained at runtime.

• •

StrictUnification: force strict unification of the gauge couplings $g_{1}=g_{2}=g_{3}$. Default false

• •

delta_mass: relative precision in the calculation of the masses. Default $10^{-6}$.

• •

n_run: maximum number of iterations for the mass calculation. Default 40.

• •

TwoLoopRGE: whether to use two-loop RGEs. Default true.

• •

Alpha: value of the fine structure constant $\alpha(0)$. Default $\tfrac{1}{137}$.

• •

gamZ: value of the decay width of the $Z$ gauge boson. Default $2.49$ GeV.

• •

gamW: value of the decay width of the $W$ gauge boson. Default $2.04$ GeV.

• •

Use_bsstep_instead_of_rkqs: use the bsstep algorithm instead of rkqs. Default false.

• •

Use_rzextr_instead_of_pzextr: use the rzextr algorithm instead of pzextr. Default false.

Note that the MSSM SPheno spectrum generator supported currently is the one immediately available after downloading SPheno, rather than including auto-genetated code from SARAH. It includes only the fixed-order Higgs mass calculation, so the hybrid EFT / fixed order calculation based on FlexibleEFTHiggs Athron:2016fuq which has been developed very recently Staub:2017jnp is not currently available in SpecBit.

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 Porod:2011nf . For any of these values, the backend convenience function run_SPheno raises an invalid_point exception, as described in gambit .

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 gambit ). 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:

⬇

• capability: unimproved_MSSM_spectrum

  function: get_CMSSM_spectrum

  options:

  mass_cut: [["h0_1",120,130],["~e-_1",100,1000]]

  mass_ratio_cut: [["~e-_1","|~chi0_1|",1,1.01]]

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 gambit , the spectrum object can then be obtained by dereferencing the safe pointer returned by the SpecBit module function, e.g.

⬇

using namespace Pipes::@\metavar{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. gambit ).

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 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 Skands:2003cj ; Allanach:2008qq 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 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 facilities, 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 ‘high-energy’.

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.

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 SpecBit-provided 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.

⬇

Spectrum spec_copy(spec); //This is a deep copy

Spectrum spec_copy = spec; //Also a deep copy

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).

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 obtain 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 spectrum-generator 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:

⬇

virtual double hard_upper() const

virtual double soft_upper() const

virtual double soft_lower() const

virtual double hard_lower() const

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++ library101010http://fthomas.github.io/slhaea/:

void writeSLHAfile(int slha_version, const str& fname)

Writes out the contents of the SubSpectrum in SLHA 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 in-memory 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 pre-existing 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 structure, 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 structure contains only two functions, getSLHAea and add_to_SLHAea. The former simply returns the contents of the structure as an SLHAea object (which can be easily written to disk, for example), while the second adds the contents of the structure 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. gambit ), 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:

⬇ // PDG code + context integer

bool has(@\metavar{tag}@, @\metavar{PDGcode}@, @\metavar{context}@)

double get(@\metavar{tag}@, @\metavar{PDGcode}@, @\metavar{context}@)

// PDG code + context integer (as std::pair)

bool has(@\metavar{tag}@, std::pair<int,int>)

double get(@\metavar{tag}@, std::pair<int,int>)

// Short name plus index (as std::pair)

bool has(@\metavar{tag}@, const std::pair<str,int>)

double get(@\metavar{tag}@, const std::pair<str,int>)

For SubSpectrum objects the above extra function signatures also exist, plus matching set and set_override functions:

⬇ // PDG code + context integer

void set(@\metavar{tag}@, @\metavar{value}@, @\metavar{PDGcode}@, @\metavar{context}@)

void set_override(@\metavar{tag}@, @\metavar{value}@, @\metavar{PDGcode}@, @\metavar{context}@)

// PDG code + context integer (as std::pair)

void set(@\metavar{tag}@, @\metavar{value}@, std::pair<int,int>)

void set_override(@\metavar{tag}@, @\metavar{value}@, std::pair<int,int>)

// Short name plus index (as std::pair)

void set(@\metavar{tag}@, @\metavar{value}@, const std::pair<str,int>)

void set_override(@\metavar{tag}@, @\metavar{value}@, std::pair<str,int>);

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).

These mechanisms together provide considerable convenience when using the Spectrum interface in module functions that do not know precisely how the backend spectrum generator has been connected to a given wrapper. As an example, any of the following calls to a Spectrum object containing a SubSpectrum that conforms to the MSSM spectrum contents (see Appendix C) will retrieve the same information (the lightest selectron/anti-selectron pole mass):

⬇

@\metavar{spectrum}@.get(Pole_Mass,"~e-",1);

@\metavar{spectrum}@.get(Pole_Mass,"~e-_1");

@\metavar{spectrum}@.get(Pole_Mass,std::make_pair("~e-",1));

@\metavar{spectrum}@.get(Pole_Mass,1000011,0);

@\metavar{spectrum}@.get(Pole_Mass,std::make_pair(1000011,0));

@\metavar{spectrum}@.get(Pole_Mass,"~e+",1);

@\metavar{spectrum}@.get(Pole_Mass,"~e+_1");

@\metavar{spectrum}@.get(Pole_Mass,std::make_pair("~e+",1));

@\metavar{spectrum}@.get(Pole_Mass,-1000011,0);

@\metavar{spectrum}@.get(Pole_Mass,std::make_pair(-1000011,0));

There is one final set of overloads to be mentioned: the safeget overloads. These act as drop-in replacements for any of the get functions, however they will throw a GAMBIT error (causing the program to halt and report a problem) if the retrieved value is not-a-number (NaN).

** 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 check-list 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 check-list 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. gambit ). These required tasks are:

Choose a SubSpectrumContents definition to which the new wrapper will conform, or write a new one (Sec. F.1). 2. 2.

Write a SpecTraits specialization to define helper types for the new SubSpectrum wrapper (Sec. F.2, but see also the extensions with the Model and Input typedefs used in Secs. F.3 and F.4). 3. 3.

Write the SubSpectrum wrapper class (Sec. F.2). This step can be further broken down into:

• •

Write the get/set interface functions to the external code, which return values conforming to the requirements of the chosen SubSpectrumContents, and connect them to the get/set interface via the fill_getter_maps and fill_setter_maps functions of the wrapper class. Sections F.2, F.3 and F.4 cover the various allowed types of interface functions and how to properly connect them to the wrapper interface.

• •

Write the special get_Model and get_Input functions in the wrapper (if needed by the chosen type of interface function; see sections F.3 and F.4)

• •

Write the RGE interface functions if RGE running capabilities are desired (Sec. F.6). 4. 4.

Write the module function(s) that construct and return a Spectrum interface object connected to the new wrapper (Sec. F.7).

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 SSDM . 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_{\text{Pl}}\approx 1.2\times 10^{19}$ GeV) has been studied extensively Sher1989 ; Elias-Miro2012 ; Alekhin2012 ; Bezrukov:2012sa ; Degrassi2012a ; Masina2013 ; Branchina2013 ; Buttazzo2013 ; DiLuzio:2014bua ; Nielsen:2014spa ; Andreassen:2014gha ; Espinosa:2015qea ; Bednyakov:2015sca (see also references within, and Refs. Lindner1986 ; Schrempp1996 ; Altarelli1994 ). Although the first bounds on the Higgs mass from vacuum stability were estimated over three decades ago Cabibbo1979 ; Hung1979 , the experimentally measured $125$ GeV Higgs Aad2012 ; Chatrchyan2012 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 Sher1989 , 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. Branchina2013 ; Branchina2015 ; DiLuzio:2015iua , using effective couplings suppressed by $1/M_{\text{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. Buttazzo:2013uya and Ref. Bednyakov:2015sca , using three-loop RGEs and two-loop threshold corrections to the quartic Higgs coupling $\lambda$ 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 Casas1995 ; Casas1996 ; Isidori2001 ; Burgess2002 ; Isidori2008 ; ArkaniHamed2008 ; Bezrukov2009 ; Hall2010 ; Ellis2009 . 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 Loebbert2015 ; Czerwinska2015 , 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 Gonderinger2010 ; Drozd2011 ; Chen2012 ; Belanger2013a ; Khan2014 ; Alanne2014 ; Han2015a ; Kanemura2015 . In these scenarios the addition 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 Khan2014 . See Ref. SSDM for a discussion of vacuum stability in a scalar singlet-extended Higgs-sector model.

Within GAMBIT we currently use FlexibleSUSY Athron:2014yba 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 Khan2014 ; Han2015a . 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 $\mu$ the Higgs potential can be approximated with an effective potential of the form

[TABLE]

for $h\gg v$, where $\mu=\mathcal{O}(h)$. 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 $\lambda(\mu)>0$ for all $\mu<M_{\text{Pl}}$ then the electroweak vacuum is the only minimum of the Higgs potential (up to $M_{\text{Pl}}$) and is therefore absolutely stable.

• •

Meta-stable: If there exists a $\mu_{0}<M_{\text{Pl}}$ such that $\lambda(\mu_{0})<0$ and the lifetime of the electroweak vacuum state is longer than the age of the universe.

• •

Unstable: If there exists a $\mu_{0}<M_{\text{Pl}}$ such that $\lambda(\mu_{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. Coleman1977 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 Kolb1990 , in which the equation of motion of the bubble of radius $R$, nucleated at time $t=0$, is $\left(|\vec{x}|^{2}-t^{2}\right)=R^{2}$.

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

[TABLE]

where $\Gamma_{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 $\phi$. For bubble nucleation in flat 4-dimensional spacetime we can take $\phi$ to be a function of $\rho=\left(\tau^{2}+|\vec{x}|^{2}\right)^{1/2}$ only (where $\tau=it$), as there exists an $O(4)$ symmetry. $B$ is given by the action Coleman1977

[TABLE]

satisfying the Euclidean equation of motion for $\phi$,

[TABLE]

The boundary conditions of this equation of motion can be combined into the single requirement that $\lim_{\rho\rightarrow\infty}\phi(\rho)=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\rightarrow-\infty$, the global vacuum at $t=0$, and bounces back to the starting point at $t\rightarrow\infty$.

The process outlined above for a simple scalar field can be applied to the Higgs field, with the use of the approximation $V(h)\approx\lambda(h)h^{4}/4$. If we perform this calculation at tree level, we can obtain an approximation for the coefficient $B$. Taking $\lambda(h)=\lambda$, $\lambda<0$, we have a solution Isidori2001 to the analogue of Eq. 4 as

[TABLE]

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

[TABLE]

The validity of the approximation $V\approx\lambda h^{4}/4$, with $\lambda<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. Isidori2001 , the bounce solution is not a constant field configuration, so requires a non-zero 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 Lee1986 .

Now that we have an expression for $B$ we need to determine the pre-factor. The explicit form of this pre-factor, $\Gamma_{0}$, was first calculated in Ref. Callan1977 by taking account of one-loop quantum corrections. Yet because the behaviour of the exponential in $\Gamma$ dominates to such a large extent, the value of $\Gamma_{0}$ need only be an approximation, so we follow the analysis of Ref. Sher1989 and determine $\Gamma_{0}$ by a dimensional reasoning.

If we take the Planck constant and the speed of light to be unity ($\hbar=c=1$), then the rate of decay per unit time per unit volume, $\Gamma/V$ has units of $[\textrm{length}]^{-4}$. Thus $\Gamma_{0}$ must have units of $[\textrm{length}]^{-4}$ or $[\textrm{energy}]^{4}$. The characteristic scale relevant in this problem is the width or height of the potential barrier, thus we take $\Gamma_{0}\approx\Lambda_{B}^{4}=1/R^{4}$, where we henceforth take $\Lambda_{B}$ to be the energy at which $\lambda(\mu)$ is at a minimum111111If the minimum value of the running quartic coupling $\lambda_{\textrm{min}}<0$ is achieved after $M_{\textrm{Pl}}$, but $\lambda(M_{\textrm{Pl}})<0$, we take $\lambda_{\textrm{min}}=\lambda(M_{\textrm{Pl}})$..

$\Gamma/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 $\Gamma/V$ by the volume of the past light cone \sim$$T_{U}^{4}, where $T_{U}$ is the age of the universe.121212This can be computed more rigorously for a standard FLRW cosmology, see Ref. Buttazzo2013 , however to the level of detail required here this result is equivalent. Thus we obtain a rate

[TABLE]

and we obtain the total expected decay rate in the past light cone

[TABLE]

Here $\lambda(\Lambda_{B})$ is the minimum value of the quartic Higgs coupling, and we have expressed the age of the Universe as $T_{u}\approx e^{140}/M_{\textrm{Pl}}$.131313Eq. 8 may be expressed in different units. For example, Ref. Sher1989 expresses the age of the universe in “units of the electroweak scale”, as $T_{u}\approx e^{101}$, where the mass of the $Z$ boson is set as $m_{z}=1$. Similarly Ref. Degrassi2014 expresses the probability in what would equivalently be called “units of the Planck scale”, but leave the original factor of $M_{\textrm{Pl}}$ in the expression; this is the style we follow.

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 $\approx 10$ billion years).141414By a decay event we mean the decay of the Universe at our position due to a decay at some point in our past light cone, and thus a decay of the observable Universe. Of course, more than one decay event does not physically make sense, so all situations where $k\geq 1$ are effectively equivalent – the Universe has decayed to the true vacuum. To model the probability that the Universe has actually decayed in the given time interval, we use a Poisson distribution, $\Theta(k;\lambda)=(\lambda^{k}/k!)e^{-\lambda}$. 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

[TABLE]

The likelihood given by Eq. 9 is typically either extremely small or exactly one, being extremely sensitive to the value of $\lambda(\Lambda_{B})$. This results in an almost step-function 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_{\mathrm{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 hard-coded to $13.6\times 10^{9}$ years. If $\lambda$ does not have a turning point before $M_{\textrm{Pl}}$, then the lifetime returned is $1\times 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_{\mathrm{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\pi)^{1/2}$ at a set of points evenly distributed in log space between the electroweak scale and the smaller of $\Lambda_{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.

⬇

- capability: print_SingletDM_spectrum

  purpose: Observable

- capability: VS_likelihood

  purpose: LogLike

- capability: expected_lifetime

  purpose: Observable

** 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. The Higgs couplings are required by HiggsBounds/HiggsSignals, which are discussed in the ColliderBit manual, Ref. ColliderBit . 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 $WW$, $ZZ$, $t\bar{t}$, $b\bar{b}$, $c\bar{c}$, $s\bar{s}$, $\tau^{+}\tau^{-}$, $\mu^{+}\mu^{-}$, $Z\gamma$ and $h\gamma$, 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).

Note that all Higgs couplings available in GAMBIT presently refer to loop-corrected values; tree-level couplings could be added if required, but should be differentiated from the present couplings by assigning them an alternative capability (e.g. Higgs_Couplings_tree_level).

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 ColliderBit 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 DarkBit 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.

Although DecayBit provides uncertainties on some decay widths and branching fractions (currently only those that we take from the PDG PDB ), these are not utilised as yet by other GAMBIT modules (indeed not even by HiggsSignals nor HiggsBounds). Future likelihood functions are expected to make use of these data, and future versions of DecayBit are anticipated to also feature calculation of uncertainties on partial and total widths of BSM states.

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.

** Supported decays**

Standard Model

DecayBit contains total widths and associated uncertainties for decays of the $W$, $Z$, $t$, $b$, $\tau$ and $\mu$, as well as their antiparticles, and the most common mesons $\pi^{0}$, $\pi^{\pm}$, $\eta$, $\rho^{0},\rho^{\pm}$ and $\omega$. It also includes partial widths to all observed, distinct final states for $W$, $Z$, $t$, $b$, $\tau$, $\mu$, $\pi^{0}$ and $\pi^{\pm}$, 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 PDB . These data are used in GAMBIT in two scenarios. The first is whenever a decay has no relevant BSM contribution, e.g. if BSM contributions are heavily loop-suppressed, or if the particle sector of the scan is simply the SM, as when scanning e.g. the DM halo model only. The second is in cases where the only significant BSM modification to the decay of the SM particle is to open a new decay channel, meaning that the new decay can simply be added to the existing SM decays. Where neither of these scenarios apply, i.e. where BSM corrections are relevant for SM $\to$ SM decays, DecayBit accounts for the corrections via an appropriately model-dependent calculation.

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. gambit .

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. YellowBook13 . Using these tables wherever possible is preferable to simply calling HDECAY, as they include additional higher-order corrections for 4-body final states obtained using PROPHECY4F Bredenstein:2006rh ; Bredenstein:2006ha .

Scalar singlet

When $m_{\scriptscriptstyle S}<m_{h}/2$, the Higgs can decay to two singlet scalars. The partial width for this process is

[TABLE]

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\rightarrow 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% Belanger:2013xza . We implement this as a complete likelihood function in DecayBit, by interpolating the function shown in Fig. 8 of Belanger:2013xza .

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. YellowBook13 . 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 Djouadi:1997yw . For Higgs masses above 1 TeV, in order to obtain finite results we disabled electroweak corrections and decays to $s\bar{s}$, $b\bar{b}$, $gg$ and $\mu^{+}\mu^{-}$ (which are negligible at large masses, where $W^{+}W^{-}$, $ZZ$ and $t\bar{t}$ final states dominate). Because of these settings, and the absence of any corrections obtained with PROPHECY4F Bredenstein:2006rh ; Bredenstein:2006ha in our supplemental calculations, small discontinuities can be expected at the transitions between different tables.

Top quark decays are available only by calling FeynHiggs. 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 as models with greater freedom rely on the more general SLHA2 format allowing full family mixing amongst the sfermions.

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, loop-induced 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 $\overline{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 $\Delta_{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

[TABLE]

we apply the correction as

[TABLE]

motivated by a similar improvement implemented in NMSDECAY NMSDECAY . Here $\Delta_{1}=\Delta_{X,1}+\Delta_{Y,1}+\ldots$ refers to the sum of 1-loop corrections to the partial widths for all 2-body final states $X$, $Y$, etc. Note that of course this expression does not actually resum higher order corrections correctly. It is simply employed as an expression that reproduces the leading order term of Eq. 11 and will be positive, so that the pathology of the negative width is avoided.

With the original definition (Eq. 11) it is thus straightforward to define a 1-loop 2-body branching fraction $\mathcal{B}$ for final state $X$ as

[TABLE]

However, this no longer follows when $\Gamma_{\rm 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 $\Delta_{1}=-1$. Whilst this is still an obvious improvement over having non-physical results for $\Delta_{1}<-1$, it is not ideal. An alternative treatment would be to use Eq. 12 whenever $\Delta_{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 $\Delta_{1}<0$. However, this would modify the standard behaviour of SUSY-HIT in the regime $-1<\Delta_{1}<0$. In the interests of keeping our modifications to external codes simple and minimal, we therefore consider it preferable to pay 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 $\overline{MS}$ 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 $\overline{MS}$ 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 significant151515We correct some minor bugs preventing compiler safety checks. 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.161616We make no attempt within GAMBIT to directly remedy any discrepancies seen between FeynHiggs and other codes Athron:2016fuq ; Staub:2017jnp ; we refer the reader to those papers, and the first GAMBIT MSSM results papers CMSSM ; MSSM for further discussion. We pass MSSM models to FeynHiggs using the function FHSetPara, directly setting the SM parameters in SLHA2 parametrisation, the $CP$-odd Higgs pole mass, the $\mu$ parameter at the SUSY scale and the $\overline{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. We set the FHSetPara arguments Qt and Qb to the renormalisation scale of the $\overline{DR}$ soft masses, so that FeynHiggs interprets stop- and sbottom-sector parameters in the $\overline{DR}$ scheme, and converts them internally to the on-shell scheme.

** 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. These are assigned names

• particle_designation_decays

and capabilities

• particle_designation_decay_rates,

where particle_designation is simply the particle name for neutral isosinglets, particle_designation = particle_name{plus*,* minus} for charged leptons and gauge bosons, particle_designation = {particle_name, particle_namebar} for quarks, and particle_designation = particle_name{0*,* plus*,* minus} for isotriplet mesons.

The function

• Ref_SM_Higgs_decays_table

constructs Reference_SM_Higgs_decay_rates from the tables of Ref. YellowBook13 (which are determined using HDECAY Djouadi:1997yw ; Spira:1997dg ; Butterworth:2010ym and PROPHECY4F Bredenstein:2006rh ; Bredenstein:2006ha ), using the pole mass provided by SpecBit or PrecisionBit (see Sec. 3.1.1). The alternative function

• Ref_SM_Higgs_decays_FH

constructs Reference_SM_Higgs_decay_rates from Higgs couplings obtained from FeynHiggs (by e.g. SpecBit::FH_Couplings; see 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\rightarrow 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.

MSSM

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.

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.171717Note that the only module functions in Table 6 with an explicit backend requirement listed are those that make use of SUSY-HIT. The corresponding FeynHiggs functions instead depend on FH_Couplings_output, which can be obtained from SpecBit’s function FH_Couplings — which exhibits the equivalent backend requirement on FeynHiggs (see Table 2). It also includes additional top decays via intermediate off-shell MSSM states, calculable via FeynHiggs. Similar to the reference functions for the SM-like Higgs described in Sec. 3.2.1, DecayBit 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 necessary 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.

** Code description and interface details**

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 PDB ), 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

// @\cpppragma{$W^+$}@ DecayTable::Entry Wplus;

// Set the total width

Wplus.width_in_GeV = 2.085;

// Set branching fraction and uncertainty

// for @\cpppragma{$W \rightarrow e\,\nu_e$}@

Wplus.set_BF(0.1071, 0.0016, "e+", "nu_e");

// Find the partial width for @\cpppragma{$W \rightarrow e\,\nu_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. gambit ). These are:

a simple argument list of strings corresponding to the final state particles’ long names (as in the example above) 2. 2.

an argument list of alternating PDG codes and context integers, 3. 3.

an argument list of alternating short name strings and integer particle indices, 4. 4.

a std::vector of long name strings, and 5. 5.

a std::vector<std::pair<int,int>> of PDG code – context 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 reminiscent 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 @\cpppragma{$W \rightarrow e\,\nu_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\rightarrow e\,\nu_{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 ColliderBit 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 (2-flavour) 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_pse-udonyms 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:

⬇ #define CAPABILITY Reference_SM_Higgs_decay_rates

START_CAPABILITY

#define FUNCTION Ref_SM_Higgs_decays_table

START_FUNCTION(DecayTable::Entry)

DEPENDENCY(mh, double)

#undef FUNCTION

#undef CAPABILITY

#define CAPABILITY Higgs_decay_rates

START_CAPABILITY

#define FUNCTION SM_Higgs_decays

START_FUNCTION(DecayTable::Entry)

DEPENDENCY(Reference_SM_Higgs_decay_rates,

DecayTable::Entry)

#undef FUNCTION

#define FUNCTION SingletDM_Higgs_decays

START_FUNCTION(DecayTable::Entry)

DEPENDENCY(Reference_SM_Higgs_decay_rates,

DecayTable::Entry)

DEPENDENCY(SingletDM_spectrum, Spectrum)

ALLOW_MODELS(SingletDM, SingletDMZ3)

#undef FUNCTION

#define FUNCTION MSSM_h0_1_decays

START_FUNCTION(DecayTable::Entry)

DEPENDENCY(SLHA_pseudonyms,

DecayBit::mass_es_pseudonyms)

BACKEND_REQ(cb_widthhl_hdec, (sh_reqd),

widthhl_hdec_type)

BACKEND_REQ(cb_wisusy_hdec, (sh_reqd),

wisusy_hdec_type)

BACKEND_REQ(cb_wisfer_hdec, (sh_reqd),

wisfer_hdec_type)

BACKEND_OPTION( (SUSY_HIT), (sh_reqd) )

ALLOW_MODELS(MSSM63atQ, MSSM63atMGUT)

#undef FUNCTION

#undef CAPABILITY

The first of these, Ref_SM_Higgs_decays_table, calculates the width and decay branching fractions of the pure SM Higgs, interpolating in the tables of Ref. YellowBook13 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\rightarrow 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 model-specific 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 gambit .

The SM reference function computes the total width and branching fractions of the Higgs as follows:

⬇ // Reference SM Higgs decays from Ref. @\cpppragma{\cite{YellowBook13}}@.

void Ref_SM_Higgs_decays_table(

   DecayTable::Entry& result)

{

using namespace

Pipes::Ref_SM_Higgs_decays_table;

// Get the Higgs pole mass

double mh = *Dep::mh;

// Invalidate point if @\cpppragma{$m_h$}@ is outside the most

// reliable part of the LHCHiggsXSWG @\cpppragma{\cite{YellowBook13}}@

// tables.

double minmass =

runOptions->getValueOrDef<double>(90.0,

"higgs_minmass");

double maxmass =

runOptions->getValueOrDef<double>(160.0,

"higgs_maxmass");

if (mh < minmass or mh > maxmass)

{

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. YellowBook13 are most reliable. The calculator and version are then set, and the total width is extracted from the Ref. YellowBook13 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:

⬇ // SM decays: Higgs

void SM_Higgs_decays (DecayTable::Entry& result)

{

namespace Pipes::SM_Higgs_decays;

result = *Dep::Reference_SM_Higgs_decay_rates;

}

The decays of the Higgs are harvested and assembled with the decays of all other particles into a single DecayTable, by the function all_decays:

⬇ #define CAPABILITY decay_rates

START_CAPABILITY

#define FUNCTION all_decays

START_FUNCTION(DecayTable)

DEPENDENCY(Higgs_decay_rates, DecayTable::Entry)

…

MODEL_CONDITIONAL_DEPENDENCY(h0_2_decay_rates,

DecayTable::Entry, MSSM63atQ, MSSM63atMGUT)

…

#undef FUNCTION

#undef CAPABILITY

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:

writing additional module functions that can generate DecayTable::Entry objects for the novel field content, 2. 2.

adding new module functions that can compute new decays of existing particles in the new model, and 3. 3.

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 gambit . 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 gambit ).

4   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 $\theta_{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 PrecisionBit’s ability to compute MSSM precision corrections and likelihoods in Sec. 4.2. Except for precision mass corrections, PrecisionBit does not (yet) support the calculation of precision corrections in non-supersymmetric models; these will be implemented as relevant when the corresponding models are implemented in GAMBIT.

Fig. 9 shows how these different components of PrecisionBit 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 $\chi^{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 (i.e. assuming that they are uncorrelated).181818Some members of the LHC Higgs Cross-section Working Group advocate a linear combination (see Sec. 12.5 of Dittmaier ), on the basis of an assumed flat prior on the theoretical systematic. However, toy Monte Carlo simulation suggests that linear combination is not a statistically valid procedure, even for this choice of prior, except in the extreme limit of complete correlation between experimental and theoretical uncertainties. Here we assume a Gaussian prior on the theoretical systematic anyway, so combination in quadrature with the (Gaussianly-distributed) experimental uncertainties is certainly well justified, as arguments associated with flat priors do not apply.

Couplings

Likelihoods associated with the Fermi coupling constant ($G_{F}$) and fine-structure constant ($\alpha_{\mathrm{em}}$) are constructed by comparing model input values with those determined from electroweak global fits, with $G_{F}=(1.1663787\pm 0.0000006)\times 10^{-5}$ GeV*-2* and $\alpha_{\mathrm{em}}(m_{Z})^{-1}=127.940\pm 0.014$ ($\overline{MS}$ scheme) PDB . Similarly, the strong coupling ($\alpha_{s}$) likelihood is calculated using a $\chi^{2}$-averaged value of $\alpha_{s}(m_{Z})=0.1185\pm 0.0005$ ($\overline{MS}$ scheme) taken from the sub-field of lattice determinations PDB . Uncertainties are taken to correspond to $1\,\sigma$ confidence intervals with no additional theoretical uncertainties. The PrecisionBit functions providing these likelihoods are summarised in Table 8.

Masses

Likelihoods associated with the $\overline{MS}$ light quark ($u,d,s$) masses at $\mu=2$ GeV, $m_{c}(m_{c})$, and $m_{b}(m_{b})$, are calculated by PrecisionBit by comparing input values from GAMBIT with the corresponding PDG values and uncertainties PDB , approximating the likelihoods as $\chi^{2}$ functions. Similarly, the likelihood corresponding to the top quark mass is evaluated using combined Tevatron and LHC measurements ATLAS:2014wva . 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 PDB . 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. PrecisionBit also provides a simple Higgs mass likelihood, using the ATLAS and CMS combination Aad:2015zhl and is used directly as the PDG average PDG15 . The module functions providing the $W$, $Z$ and $h$ mass likelihoods are summarised in Table 8. Unlike the other SM particle masses, the model-dependent $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 $\chi^{2}$ approximation to the mass likelihoods. By summing in quadrature, we assume that the experimental and theoretical errors are uncorrelated, which is reasonable given that they come from entirely different sources. The more detailed Higgs likelihoods of ColliderBit ColliderBit 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. SSDM .

** 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 Heinemeyer:1998yj , SuperIso Mahmoudi:2008tp and GM2Calc gm2calc 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 FlavBit .

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_{\mu}$ to $(g-2)_{\mu}$. 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, $\alpha$, in the Thompson limit and another $\alpha$ that includes SM fermion contributions to the on-shell photon vacuum polarization. These definitions differ from the $\overline{MS}$ $\alpha(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.

PrecisionBit also contains an experimental interface to the proprietary code SUSY-POPE Heinemeyer:2006px ; Heinemeyer:2007bw , but this is not recommended for general use. In addition, a new FeynHiggs routine FHEWPO, based on work used for Refs. Heinemeyer:2013dia ; Stal:2015zca , has very recently been released for beta testing in FeynHiggs 2.13.0-beta. We understand that this beta is now to be considered production-stable,191919T. Hahn, personal communication and so expect to support it in a future version of GAMBIT. Further calculations of the EWPO may in future become available through the FlexibleTools collaboration behind FlexibleSUSY or through use of auto-generation software, which is now being developed with GAMBIT.

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

Electroweak precision observables are well measured observables from the electroweak sector of the SM. Due to the precise measurements they can provide constraints on BSM physics. For a pedagogical introduction see Ref. Matchev:2004yw , where $Z$-pole and non $Z$-pole sets of electroweak precision observables are listed in Tables 1 and 2, respectively. 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}\theta_{W,\mathrm{eff}}$, the departure from 1 of the ratio of the Fermi constants implied by neutral and charged currents $\Delta\rho$, and the muon anomalous magnetic moment $a_{\mu}$ (also referred to as $g-2\equiv 2a_{\mu}$). For each of these observables, PrecisionBit calculates likelihoods by combining experimental and theoretical uncertainties in quadrature (i.e. assuming that theoretical and experimental errors are not correlated), 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.

PrecisionBit also offers two module functions able to provide the predicted SM contribution to $g-2$. These are listed along with the likelihood functions in Table 11. The first is calibrated to $e^{+}e^{-}$ data, giving $a_{\mu,\mathrm{SM}}=(11659180.2\pm 4.9)\times 10^{-10}$ 1010.4180 , whereas the second comes from $\tau^{+}\tau^{-}$ data, and corresponds to $a_{\mu,\mathrm{SM}}=(11659189.4\pm 5.4)\times 10^{-10}$ 1010.4180 .

For the $g-2$ likelihood, we combine errors from the predicted SM contribution in quadrature with the theoretical error arising from the BSM contribution, as these two theory errors are entirely independent of each other, and therefore uncorrelated. In the likelihood function, we compare the sum of the SM and BSM contributions to the experimental measurement $a_{\mu}=(11659208.9\pm 6.3)\times 10^{-10}$ PDG10 ; gm2exp , where the error is the sum in quadrature of statistical ($5.4\times 10^{-10}$) and systematic ($3.3\times 10^{-10}$) uncertainties.

Theoretical calculations of the MSSM contribution to $a_{\mu}$ 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 its prediction, using the magnitudes of the two-loop Barr-Zee corrections barrzee to estimate the magnitude of neglected higher-order contributions gm2calc ; PrecisionBit adopts this estimate directly when employing GM2Calc. For FeynHiggs and SuperIso, no such estimate is available. Following the discussion in Ref. stockinger06 , we assign a theoretical uncertainty of either 30% or $6\times 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}\theta_{W,\mathrm{eff}}$ and $\Delta\rho$ in the MSSM are currently provided exclusively by FeynHiggs. We assume a theoretical error, $\sigma^{\mathrm{theo}}_{\sin^{2}\theta_{W,\mathrm{eff}}}$, of $12\times 10^{-5}$ on $\sin^{2}\theta_{W,\mathrm{eff}}$ hep-ph/0412214 . The fractional MSSM corrections to both $\sin^{2}\theta_{W,\mathrm{eff}}$ and $m_{W}$ are approximately proportional to $\Delta\rho$ at tree level, with a constant of proportionality quite close to one in magnitude (within a few tens of percent; hep-ph/0412214 ). Combined with the fact that $\Delta\rho\ll 1$ in general, to a good approximation the theoretical error on its value can be estimated from the fractional uncertainty on $\sin^{2}\theta_{W,\mathrm{eff}}$ or $m_{W}$; we therefore adopt a theory error on $\Delta\rho$ of

[TABLE]

where $\sigma^{\mathrm{theo}}_{m_{W}}=10$ MeV hep-ph/0412214 .

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:

• make_MSSM_precision_spectrum_H_W(),
• make_MSSM_precision_spectrum_W(), and
• make_MSSM_precision_spectrum_none().

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 Higgs-sector 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 FeynHiggs (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 hep-ph/0412214 (as discussed in the context of the theory error on $\Delta\rho$ in Sec. 4.2.2).202020It is however worth noting that this is in fact no smaller than the value that we assign to the uncertainty on the $W$ mass calculation performed by FlexibleSUSY, so the replacement has no effect if FlexibleSUSY has been used to generate the original unimproved_MSSM_spectrum.

For the MSSM Higgs bosons ($A^{0}$, $h^{0}$, $H^{0}$ and $H^{\pm}$), 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.

• 2:

Keep the masses unchanged from the original input spectrum.

• 3:

Assign the mean of the masses found in the original spectrum and the ones provided by the precision calculator.

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 ($\Delta_{s}$), the error determined by the precision calculator ($\Delta_{p}$), and the difference of central values from the two calculations ($\Delta_{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 $\sigma_{+}$ 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 $\sigma_{-}$ 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: $\sigma_{+}^{2}=\sigma_{-}^{2}=\Delta_{s}^{2}+\Delta_{p}^{2}+\Delta_{g}^{2}$.

• 2:

RACC, with $\Delta_{s}$ added to the error associated with the distance to the spectrum generator value (the spectrum-generator ‘edge’), and $\Delta_{p}$ added to the precision-calculator edge. This is the default.

• 3:

RACC, with $\Delta_{g}/2$ added to both $\sigma_{+}$ and $\sigma_{-}$.

• 4:

RACC, with $\Delta_{g}/2$ added at the spectrum-generator edge, and $\Delta_{p}$ added at the precision-calculator edge.

• 5:

RACC, with $\Delta_{g}/2$ added at the precision-calculator edge, and $\Delta_{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.

5   Examples

** 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 gambit . DecayBit_SingletDM does the same for ten randomly-selected 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 parametrisation. 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 3bithit.out.slha.

The basic spectrum generation in 3-BIT-HIT is done using FlexibleSUSY via SpecBit, with the Higgs and $W$ masses replaced by precision calculations performed with FeynHiggs by PrecisionBit. The resulting SM-like Higgs mass and uncertainties emitted in the final SLHA file are based on the default settings in PrecisionBit::make_MSSM_precision_spectrum_H_W():

⬇

Higgs_predictions_source = 1

Higgs_predictions_error_method = 2

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

6   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.

Acknowledgements.

We thank all our colleagues within GAMBIT for helpful discussions, especially Nazilla Mahmoudi for her comments on the muon $g-2$ functionality of SuperIso. We would also like to thank Margarete Mühlleitner and Ramona Gröber for responding to questions regarding SDECAY, and Sven Heynemeyer and Sebastian Paßehr for doing so on FeynHiggs. 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. GAMBIT has been supported by STFC (UK; ST/K00414X/1, ST/P000762/1), the Royal Society (UK; UF110191), Glasgow University (UK; Leadership Fellowship), the Research Council of Norway (FRIPRO 230546/F20), NOTUR (Norway; NN9284K), the Knut and Alice Wallenberg Foundation (Sweden; Wallenberg Academy Fellowship), the Swedish Research Council (621-2014-5772), the Australian Research Council (CE110001004, FT130100018, FT140100244, FT160100274), The University of Sydney (Australia; IRCA-G162448), PLGrid Infrastructure (Poland), Polish National Science Center (Sonata UMO-2015/17/D/ST2/03532), the Swiss National Science Foundation (PP00P2-144674), the European Commission Horizon 2020 Marie Skłodowska-Curie actions (H2020-MSCA-RISE-2015-691164), the ERA-CAN+ Twinning Program (EU & Canada), the Netherlands Organisation for Scientific Research (NWO-Vidi 680-47-532), the National Science Foundation (USA; DGE-1339067), the FRQNT (Québec) and NSERC/The Canadian Tri-Agencies Research Councils (BPDF-424460-2012).

Appendices

Appendix 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 code-focussed description in the main GAMBIT paper gambit . 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. DarkBit for an example in GAMBIT). We therefore begin with the SM itself.

Standard Model

The SM describes the interactions of all observed fundamental particles. These interactions are invariant under the symmetries of the gauge group,

[TABLE]

The Lagrangian density is given by,

[TABLE]

where $u_{Ri}$, $d_{Ri}$ and $e_{Ri}$ represent the right-handed up-type quark, down-type quark and lepton. Similarly,

[TABLE]

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)_{\mathrm{L}}$ indices, and $\epsilon_{ab}$ is the anti-symmetric tensor.

The $B_{\mu\nu}$, $W^{A}_{\mu\nu}$ and $G^{B}_{\mu\nu}$ are field strength tensors for the $B_{\mu}$, $W^{A}_{\mu}$ and $G^{B}_{\mu}$ fields respectively, which are respectively associated with the $U(1)_{\mathrm{Y}}$, $SU(2)_{\mathrm{L}}$ and $SU(3)_{\mathrm{C}}$ gauge groups. The kinetic terms for the matter fields include $\not{D}=D^{\mu}\gamma_{\mu}$, where $D_{\mu}$ is the usual covariant derivative, which, for example, takes the following form for all fields transforming non-trivially under both $SU(2)_{\mathrm{L}}$ and $SU(3)_{\mathrm{C}}$:

[TABLE]

Here $g_{1}$ and $Q_{Y}$ are the GUT-normalised gauge coupling and GUT-normalised charges212121Hypercharge is often written without GUT normalisation, such that the charge for a field $\phi$ is $Y^{\phi}=\sqrt{5/3}Q_{Y}^{\phi}$ and will be written in the covariant derivative with the gauge coupling $g^{\prime}=\sqrt{3/5}g_{1}$. of $U(1)_{Y}$. Similarly $g_{2}$ and $T^{A}$ are the $SU(2)_{\mathrm{L}}$ gauge couplings and generators respectively, and $g_{3}$ and $\lambda^{B}$ are the $SU(3)_{\mathrm{C}}$ gauge couplings and generators respectively.

Electro-weak symmetry is broken when the Higgs field develops a VEV, $v$,

[TABLE]

so that field may be rewritten as

[TABLE]

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.

The SM VEV gives masses to the SM particles. The $Z$ and $W$ gauge bosons obtain tree-level masses

[TABLE]

The weak mixing angle $\theta_{W}$ is given by

[TABLE]

though it is important to note that the last equality depends on the tree-level relations given in Eq. 21, 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,

[TABLE]

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.

Minimal Supersymmetric Standard Model

The MSSM superpotential is given by,

[TABLE]

where a hat indicates a superfield and $i$, $j$ specify the generation, while $a$, $b$ are $SU(2)_{\mathrm{L}}$ indices. This includes the Higgsino mixing parameter $\mu$, which appears as a dimension one coefficient of the up- and down-type Higgs superfields,

[TABLE]

The chiral superfields have dimensionless up- and down-type Yukawa couplings, $\mathbf{Y}_{u}$ and $\mathbf{Y}_{d}$ with the Higgs superfields and the chiral superfields for the right-handed quarks. Specifically left-handed quarks, $\hat{Q}$, interact with both the chiral superfields for the right-handed up quarks, $\hat{U}^{c}$, and the up-type Higgs superfield through $\mathbf{Y}_{u}$ and with down quarks, $\hat{D}^{c}$, and the down-type Higgs superfields through $\mathbf{Y}_{d}$. Similarly the chiral superfield for left-handed leptons $\hat{L}$ has Yukawa interactions, $\mathbf{Y}_{e}$ with the chiral superfields for the right-handed leptons, $\hat{E}^{c}$ and $\hat{H}_{d}$.

Terms which break supersymmetry without introducing quadratic divergences are given in the soft SUSY breaking Lagrangian,

[TABLE]

As in the previous equation we have indices $i$, $j$ for generation and $a$, $b$ for $SU(2)_{\mathrm{L}}$ and in addition $A\in\{1,2,3\}$ is the index for $SU(2)_{\mathrm{L}}$ gauge generators and $B\in\{1,2,3...8\}$ is the index for the $SU(3)_{\mathrm{C}}$. The superpartners of the $U(1)_{Y}$ gauge field, $\tilde{B}^{0}$, the $SU(2)_{\mathrm{L}}$ gauge fields, $\tilde{W}_{A}$, and the $SU(3)_{\mathrm{C}}$ gauge fields, $\tilde{g}_{B}$ get soft masses, $M_{1}$, $M_{2}$ and $M_{3}$ respectively, as well as potential contributions from $M_{1}^{\prime}$, $M_{2}^{\prime}$ and $M_{3}^{\prime}$. The Higgs scalar fields, $H_{u}$ and $H_{d}$ receive soft scalar mass contributions, $m^{2}_{H_{u}}$, $m^{2}_{H_{d}}$ respectively and are mixed by the soft bilinear mass squared, $b$. The left-handed squarks, $\tilde{Q}_{i}$ and sleptons, $\tilde{L}_{i}$ and the right handed up-type squarks, $\tilde{u}_{\mathrm{Ri}}$, down-type squarks, $\tilde{d}_{\mathrm{Ri}}$, and charged sleptons $\tilde{e}_{\mathrm{Ri}}$ all have soft scalar masses, $(\mathbf{m}^{\mathbf{2}}_{Q})_{ij}$, $(\mathbf{m}^{\mathbf{2}}_{L})_{ij}$, $(\mathbf{m}^{\mathbf{2}}_{u})_{ij}$, $(\mathbf{m}^{\mathbf{2}}_{d})_{ij}$ and $(\mathbf{m}^{\mathbf{2}}_{e})_{ij}$ respectively. Finally we have soft trilinears for the up-type squarks and Higgs, $(\mathbf{T}_{u})_{ij}$, down-type squarks and Higgs, $(\mathbf{T}_{d})_{ij}$, and the sleptons and down-type Higgs, $(\mathbf{T}_{e})_{ij}$ and also the terms that are are soft when there is no singlet in the model, $(\mathbf{C}_{u})_{ij}$, $(\mathbf{C}_{d})_{ij}$ and $(\mathbf{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,

[TABLE]

In the MSSM the neutral components of the up- and down-type scalar Higgs fields, $H_{u}^{0}$ and $H_{d}^{0}$ respectively, develop VEVs,

[TABLE]

and the bilinear soft term, $b$, mixes them, so that they must be rotated to obtain the physical mass eigenstates, $h_{i}$

[TABLE]

where the mixing matrix is commonly expressed in terms of a mixing angle $\alpha$,

[TABLE]

Similarly the pseudoscalar Higgs state, $A$, is formed from the imaginary parts of the Higgs scalar fields

[TABLE]

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,

[TABLE]

giving a massless Goldstone boson and an expression for the pseudoscalar mass,

[TABLE]

The charged Higgs mass eigenstate, $H^{+}$ and charged Goldstone boson, $G^{+}$ are related to the charged scalar Higgs fields through,

[TABLE]

The squarks are rotated into the mass eigenstates as,

[TABLE]

and the sleptons and sneutrinos are rotated as,

[TABLE]

The neutralinos are mass eigenstates from combinations of the Bino $\tilde{B}^{0}$, neutral Wino $\tilde{W_{3}}$ and the up- and down-type Higgsinos, $\tilde{H}^{0}_{u}$ and $\tilde{H}^{0}_{d}$,

[TABLE]

where we follow the SLHA Skands:2003cj convention and require the mixing matrix, $Z_{N}$, to be real, which means the masses of the neutralinos may be negative. Note that this is the form of the masses and mixings that can be accessed from the spectrum object in GAMBIT and used for passing into interfaces to various backend codes that need the spectrum as an input. We follow the SLHA convention of negative masses because this ensures the maximum number of codes can be supported without conversion during interfacing. For any codes or internal calculations which expect neutralino masses and mixings in another convention, this must be converted in the interface.

The mass matrix of the charginos is diagonalised through a biunitary transformation, such that the two component chargino mass eigenstates are given by,

[TABLE]

[TABLE]

Scalar Singlet Dark Matter Model

The scalar singlet model is a minimalistic extension of the SM, with one new scalar charged only under a $\mathbb{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,

[TABLE]

where $\mu_{S}^{2}$ is a mass squared term for the singlet, $\lambda_{h\scriptscriptstyle S}$ the Higgs-portal coupling and $\lambda_{\scriptscriptstyle S}$ the quartic self coupling.

The singlet mass has an additional contribution from EWSB through the Higgs-portal coupling, giving a mass,

[TABLE]

where $v$ is the SM Higgs VEV in Eq. 19. As the $\mathbb{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 $\lambda_{S}>0$.

Appendix 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 directly 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 gambit .

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 FlexibleSUSY Athron:2014yba internally implemented as the de facto spectrum generator for MSSM models, but also allows the use of the backended spectrum generator SPheno Porod:2003um ; Porod:2011nf . 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 Allanach:2001kg ; Allanach:2009bv ; Allanach:2011de ; Allanach:2014nba or SuSpect Djouadi:2002ze . 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 gambit ).

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 gambit , 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 easily 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 gambit , 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. gambit , 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 output 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.

Input and output, warnings and errors

SPheno follows the SLHA conventions Skands:2003cj ; Allanach:2008qq 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 parameters are initialised in the convenience function InitializeStandardModel using the structure 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 TEIN 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 post-generation 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 Porod:2011nf .

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:

⬇

#define FUNCTION get_MSSM_spectrum_SPheno

START_FUNCTION(Spectrum)

ALLOW_MODELS(CMSSM, MSSM63atMGUT, MSSM63atQ)

DEPENDENCY(SMINPUTS, SMInputs)

BACKEND_REQ(SPheno_MSSMspectrum,

(libSPheno), int, (Spectrum&, const Finputs&) )

BACKEND_OPTION((SPheno, 3.3.8), (libSPheno))

#undef FUNCTION

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 FlavBit , DecayBit and ColliderBit ColliderBit .

Appendix 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,

⬇

@\metavar{subspectrum}@.get(@\metavar{tag}@,@\metavar{label}@)

@\metavar{subspectrum}@.get(@\metavar{tag}@,@\metavar{label}@,@\metavar{index 1}@)

@\metavar{subspectrum}@.get(@\metavar{tag}@,@\metavar{label}@,@\metavar{index 1}@,@\metavar{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 [math], $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 GAMBIT 1.0.0, and describe each of them in detail in the table referenced beside each list item.

Appendix 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 of SpecBit. Wrappers based on simple parameter containers with no RGE facilities are listed in Table D, 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.