anyH3: precise predictions for the trilinear Higgs coupling in the Standard Model and beyond

The trilinear Higgs coupling $\lambda_{hhh}$ of the detected Higgs boson is an important probe for physics beyond the Standard Model. Correspondingly, improving the precision of the theoretical predictions for this coupling as well as the experimental constraints on it are among the main goals of particle physics in the near future. In this article, we present the public $\mathtt{Python}$ code $\mathtt{anyH3}$, which provides precise theoretical predictions for $\lambda_{hhh}$. The program can easily be used for any renormalisable model, where for the input the $\mathtt{UFO}$ format is adopted. It allows including corrections up to the full one-loop level with arbitrary values of the external squared momenta and features a semi-automatic and highly flexible renormalisation procedure. The code is validated against known results in the literature. Moreover, we present new results for $\lambda_{hhh}$ in models which so far have not been investigated in the literature.


Introduction
The discovery of a Higgs boson at the CERN LHC [1,2] has confirmed that the Higgs potential plays a crucial role in the electroweak symmetry breaking (EWSB). While the measured properties of this Higgs boson are so far compatible with the predictions of the Standard Model (SM) within the experimental and theoretical uncertainties, the structure of the Higgs sector and of its potential remain to be determined. Furthermore, in spite of the successes of the SM, it is clear that new, Beyond-the-Standard-Model (BSM) physics is needed to address deficiencies of the SM -such as for instance the lack of an explanation for the baryon asymmetry of the Universe.
In this context, a key quantity to investigate is the trilinear Higgs coupling λ hhh . This coupling determines the shape of the Higgs potential away from the electroweak (EW) minimum and in turn controls the nature and strength of the EW phase transition (EWPT). For instance, a strong first-order EWPT, which is a requirement [3] for the scenario of EW baryogenesis [4,5], is typically associated with a sizeable deviation of λ hhh from its SM prediction [6,7]. Even beyond its crucial role in the context of the EWPT, λ hhh provides a unique opportunity to find signs of BSM physics arising from extended Higgs sectors. In particular, the loop contributions in models with additional scalars can cause the trilinear Higgs coupling to deviate by up to several hundred percent from its prediction in the SM if there is a substantial mass splitting between the BSM mass scales. This was found first at the one-loop level for the case of a Two-Higgs-Doublet Model (THDM) in Refs. [8,9] but is now known to occur for a wide range of BSM models. The genuine physical nature of these effects was confirmed in Refs. [10,11], where next-to-leading order (NLO), i.e. two-loop, corrections to those large one-loop effects were investigated and found to obey the expected perturbative behaviour. Unlike what is the case for most couplings of the SM-like Higgs boson at 125 GeV, large deviations in λ hhh are possible even in scenarios where all its couplings are very close to the SM values at tree level, such as in aligned scenarios [12]. Meanwhile, it has recently been shown in Ref. [13] that the experimental limits (discussed in further detail below) have become sufficiently strong to probe these potentially large loop effects, and thus the comparison of the predictions for λ hhh with the latest experimental bounds constitutes a powerful new method for constraining the parameter space of BSM theories. In this context it is important to keep in mind that λ hhh cannot be directly measured experimentally. The crucial experimental quantity where λ hhh enters at leading order is the process of Higgs pair production. The computation of the trilinear Higgs coupling λ hhh constitutes a necessary intermediate result for the prediction of di-Higgs boson production. In fact, for the case where λ hhh receives large loop corrections, the additional contributions to the di-Higgs production process may be of sub-leading order [13]. More generally, computations of trilinear Higgs couplings are also important for investigating decays of BSM Higgs bosons and BSM decays of the SM-like Higgs boson.
The present experimental information on the trilinear Higgs coupling λ hhh is by far not as precise as what has been achieved for other couplings of the Higgs boson [14,15]. Indeed, the current best limits on λ hhh were obtained by the ATLAS collaboration using a combination of data from searches for (non-resonant) di-Higgs production and from the experimental results for single-Higgs production processes; they bound the ratio κ λ , defined as to be within the range −0.4 < κ λ < 6.3 at the 95% confidence level (C.L.) [16,17]. In Eq. (1), (λ SM hhh ) (0) denotes the tree-level prediction for the trilinear coupling in the SM. The CMS collaboration has obtained similar results [15], namely −1.24 < κ λ < 6.49. The quoted limits were obtained under the assumption that besides a variation of κ λ all other couplings entering the analyses are fixed to their SM values. The current experimental limits leave ample room for BSM deviations, which would so far remain unobserved, but could be accessed in the foreseeable future, given the expected prospects for probing λ hhh at the LHC and future colliders -see Ref. [18] for a review. Specifically, at the high-luminosity upgrade of the LHC (the HL-LHC), the projection for κ λ at 95% C.L. is 0.1 < κ λ < 2.3 [19]. At a future e + e − linear collider and a 100-TeV hadron collider it is expected that κ λ can be determined at the level of O(10%) [18,[20][21][22][23][24]. It should be noted that these projections were obtained under the assumption that κ λ = 1 is realised in nature and may significantly change if the actual value of κ λ is different. In particular, for an enhanced value of κ λ the prospects for extracting λ hhh from the process e + e − → Zhh at a linear collider with about 500 GeV would improve, while as a consequence of destructive interference contributions the prospects at the HL-LHC would deteriorate, see e.g. Refs. [25][26][27].
As stressed above, the most direct probe of the trilinear Higgs coupling are searches for di-Higgs production, because this process involves λ hhh already at the leading order (LO). Single-Higgs production involves contributions of the trilinear Higgs coupling starting at the next-to-leading order (NLO) and EW precision observables at the next-to-next-to-leading order (NNLO) -see for instance Refs. [28,29]. Of course, a general analysis should not be restricted to the case where BSM contributions enter exclusively via the trilinear Higgs coupling. On the other hand, in scenarios where large loop corrections to λ hhh constitute the leading contributions to di-Higgs production, an effective coupling approach where the dominant corrections are incorporated into κ λ provides a convenient framework to efficiently constrain BSM models with available experimental results, as discussed e.g. in Ref. [13]. We will discuss in this paper in more detail the applicability of experimental constraints set on λ hhh .
In this work, we present the Python package anyH3, which takes a big step forward in facilitating the prediction of λ hhh . 1 anyH3 allows the analytic and/or numerical computation of the trilinear Higgs coupling for general renormalisable theories to full one-loop order. For user convenience, the model definitions needed in anyH3 in order to enable the application of generic results to specific theories can be provided in the form of the widely-employed UFO format [57,58]. anyH3 also offers a high level of flexibility in the renormalisation schemes used in calculations -with pre-defined commands for standard scheme choices and the additional possibility for the user to define other choices of counterterms. Furthermore, the tool allows the user to modify the treatment of tadpole contributions (for recent discussions see e.g. Refs. [59][60][61][62]). anyH3 is part of the wider anyBSM framework, where developments for further (pseudo-)observables are foreseen in the future. An additional anyBSM feature that is already available is the module anyPerturbativeUnitarity, which allows efficient and reliable verifications of perturbative unitarity constraints (at leading order, and in the high-energy limit).
This paper is organised as follows: we start by discussing in Section 2 the main elements of our automated computation of λ hhh , as well as the interpretation of the obtained results. Next, we present in Section 3 the workflow of anyH3 before presenting a brief tutorial of the program in Section 4. In Section 5, we discuss the cross-checks that were performed for the various models that are installed along with anyH3. Finally, we present in Section 6 examples of applications of anyH3 with an emphasis on new results. We summarise our results in Section 7. A number of Appendices provide additional details on the program and the considered models: Appendix A contains our conventions for general renormalisable models and for the generic expressions included in anyH3 as well as the conventions or restrictions on UFO model files; Appendix B discusses the various ways of generating compatible UFO models for new BSM theories; Appendix C presents the different models discussed in this paper, including details on renormalisation prescriptions and the treatment of tadpole contributions; Appendix D lists modifications of the UFO standard that are applied by anyBSM internally; Appendix E explains the caching that is available in anyH3, while Appendix F describes the Python interface pyCollier to the Fortran library COLLIER employed for computing loop functions numerically.

A generic approach to the trilinear Higgs coupling
In the following we outline the steps for the calculation of the trilinear Higgs coupling where the SM-like Higgs boson h appears at each of the three external legs. A generalisation to trilinear Higgs couplings involving one or more BSM Higgs boson(s) is left for future work.
All calculations performed in the program are based on results for general renormalisable theories. These generic results were obtained by the following steps: • For the different types of contributions entering the calculation, as specified in Eq. (2) below, all possible Feynman diagram topologies were identified.
• For each of these topologies, all possible insertions of generic fields (scalars, fermions, vector bosons, ghosts) were processed.
The resulting generic expressions were hard-coded into the Python program code to be applied to a specific model upon run-time. The steps described above include the introduction of several conventions in how the generic Lagrangian and its resulting Feynman rules are written. For instance, all fermionfermion-scalar operators, F 1 F 2 S 3 , are written in terms of left-and right-handed projectors, (c 123 L P L + c 123 R P R )F 1F2 S 3 , where the explicit form of c 123 L,R (depending on the corresponding operator within the considered model) is yet unspecified. A detailed description of all used conventions is given in Appendix A. The generic results, which are expressed in terms of generic couplings, have been implemented into the Python code. Upon run-time of the program, the couplings of the specified model are mapped onto the generic Lagrangian allowing one to directly obtain results for all contributing Feynman diagrams. A similar approach is followed for example in SARAH [67][68][69][70][71] and TLDR [72].
Currently anyBSM contains generic results for the scalar three-point function (where in the present implementation the three external scalars are assumed to be identical), the scalar two-point function, the scalar one-point function, the vector-boson two-point function, and the vector-boson-scalar two-point function. An overview of all topologies can be found in Appendix A.2. Making use of these building blocks, the trilinear Higgs coupling λ hhh is calculated by the sub-module anyH3 at the one-loop level, λ hhh = −Γ hhh (p 2 1 , p 2 2 , p 2 3 ) = λ (0) hhh + δ (1) genuine λ hhh + δ (1) tadpoles λ hhh + δ (1) WFR λ hhh + δ (1) CT λ hhh , (2) whereΓ hhh is the renormalised Higgs-boson three-point function. The superscripts indicate the loop order. The last term of Eq. (2) denotes the counterterm contribution (see below), while the second-last term encodes the contribution from external leg corrections. As indicated in Eq. (2), anyH3 is able to handle arbitrary external momenta. For all the one-loop pieces appearing in Eq. (2), we work in anyH3 only with UV-finite parts, unless otherwise specified. Explicit checks of UV-finiteness can be performed in anyH3, but have been done separately.
We also note that anyH3 allows independent calculations of self-energies and tadpoles. The necessary generic results were derived following the same steps as described above.

Renormalisation
Renormalisation is a fundamental ingredient of loop calculations. Minimal subtraction schemes like MS are arguably the easiest schemes to automate since their implementation basically boils down to setting the divergent parts of the appearing loop integrals to zero.
However, for many processes it is known that MS schemes may result in undesirable features like artificially large loop corrections or gauge dependencies. Also in the specific context of the trilinear Higgs coupling it has been observed that a potentially large part of the one-loop corrections can be absorbed into the Higgs-boson mass, which appears at the tree level, for instance by renormalising it in the on-shell (OS) scheme [11,31,32].
For this reason, anyH3 allows the specification of different renormalisation schemes. All bosonic (i.e., of particles with spin zero or spin one) masses appearing in the tree-level expression for the trilinear Higgs coupling can optionally be renormalised in the OS scheme. Moreover, also the vacuum expectation value (VEV) entering at lowest order can be renormalised in the OS scheme (see Appendix C.1.2 for more details). In addition, anyH3 allows the user to define custom counterterms, which are then included in the calculation of the trilinear Higgs coupling. In summary, the counterterm contribution to λ hhh reads where h is used to denote the SM-like Higgs boson, m i denotes all scalar or vector boson masses appearing in λ (0) hhh , and v is used to denote the electroweak VEV. The notations δ CT x denote the one-loop counterterms for the parameters x. If one of the masses m i is chosen to be renormalised in the OS scheme, its counterterm is determined via where Σ ii is the self-energy of the scalar/vector particle i (with mass m i ), which is defined according to the conventions of e.g. Refs. [72,73]. 2 For the determination of the electroweak VEV counterterm, we refer to Appendix C.1.2.

Tadpole contributions
Besides the renormalisation of the masses and the electroweak VEV appearing at the tree level, also tadpole contributions need to be taken into account. Since the UFO standard does not provide a unified notation to store information about the minimisation of the Higgs potential 3 , we use as default setting the Fleischer-Jegerlehner tadpole scheme of Ref. [74] (using MS tadpole counterterms). As a consequence, all tadpole diagrams have to be calculated explicitly -see Appendix C.1.1 for a detailed discussion. This does not only include explicit tadpole contributions to λ hhh , denoted as δ (1) tadpoles λ hhh , but also to δ CT λ hhh and δ WFR λ hhh . It should be stressed that the treatment of the tadpole contributions can be adapted by the user. In particular, it is possible to avoid the explicit appearance of tadpole diagrams by an appropriate choice of δ (1) custom-CT λ hhh using e.g. OS tadpole counterterms. 4 2 We note that this is also the sign convention employed by anyH3 internally. For separate calculations of self-energies a flag can be used to switch the overall sign convention for the results, see the online documentation for more details. 3 The dependence of masses and vertices on tadpole terms is at present not stored in the UFO model files. 4 Note that δ (1) custom-CT λ hhh can be defined in terms of one-, two-, or three-point functions (and derivatives thereof), which are computed automatically by the code.

External leg corrections
anyH3 also includes external leg corrections to ensure the proper normalisation of the external scalars. These corrections are given by where the prime indicates a derivative with respect to the external momentum squared. The second line of this equation serves to define the notations δ (1) Z(p 2 ), which we will employ later in this paper. For the on-shell case, p 2 i = m 2 h , the first term on the right-hand side yields the usual LSZ factor as it occurs for the case without mixing between different Higgs bosons, while the second term accounts for the contributions from possible scalar mixing effects on the external legs. As explained above, the self-energies appearing in Eq. (5) are meant to contain only the UV-finite contributions. It should be noted that anyH3 by default evaluates the external leg corrections at the same momenta as the vertex corrections. In order to implement different choices of the field renormalisations together with their appropriate wave function normalisation contributions, one can alternatively choose to turn off the automatic calculation of external-leg corrections and re-introduce the corresponding contributions in δ (1) custom-CT λ hhh . This is in particular needed for the case where the result for the trilinear Higgs coupling obtained with anyH3 is meant to be incorporated into the prediction for the cross section for di-Higgs production. In the di-Higgs production process the trilinear Higgs coupling enters with two on-shell external legs, while the third leg is an off-shell internal line of the amplitude for di-Higgs production, see Section 6.5 and Appendix C.2 below for more details.

Interpretation of the result for λ hhh
The trilinear Higgs coupling λ hhh itself is not a physical observable. Its experimental determination from a physical process (typically di-Higgs production, but via higher-order contributions also single-Higgs production provides some sensitivity) relies on assumptions on the other couplings and particles involved in the processes. The current limits on λ hhh from ATLAS [16] and CMS [15] were obtained under the assumption that all other Higgs couplings that are relevant for the respective processes have the SM values and that no other BSM particles contribute to these processes. Moreover, λ hhh has been treated as a constant which does not depend on the inflowing momenta. When comparing the predictions obtained with anyH3 in the considered model to the experimental limits, the user of anyH3 should ensure that these assumptions are fulfilled sufficiently well.
Alternatively, the output of anyH3 can be used as input for other codes calculating di-Higgs production cross sections. In this case, the user can choose to include the momentum dependence and switch-off the external-leg corrections for the internal Higgs propagator (an explicit example for this case is discussed in Fig. 13 and also in Appendix C.2).

User-and Program-flow
The main objective of this work is to provide one-loop corrections to the trilinear Higgs coupling in wide classes of BSM models. A typical work-flow of how this is organised is shown in Fig. 1. We distinguish two major sections for demonstrative purposes: A) User input and B) the actual program flow. User input is required in section A). In addition, it is also possible to control each of the steps discussed in B) by using the anyBSM library. The latter will be discussed in more detail below. At the user level (A) several inputs are required: • A.1) The model-specific information such as particle content and Feynman rules relies on the UFO standard. In order to obtain a UFO description of the model of interest one can use SARAH, FeynRules, or UFO files obtained with any other tool. For the latter two cases, we provide a converter, which is discussed in Appendix A.1, while in the case of SARAH the conventions for all relevant Lorentz structures match the anyBSM conventions. A detailed description of the UFO format can be found in Refs. [57,58].
• A.2) In order to renormalise the relevant input parameters consistently, the tool needs to know which of the particles defined in the UFO model correspond to the SM particles. This information is specified in an auxiliary file called schemes.yml.
• A.3) Once the SM parameters and particles (along with their masses) are identified, one needs to specify in which renormalisation scheme they are given. This is also done in the file schemes.yml. An example specification of this file is shown in Section 4.4.
• A.4) Numerical values for all input parameters. Additionally, the UFO model may provide analytic relations between input parameters (so called "external" parameters) and e.g. Lagrangian parameters or mixing angles (so called "internal" parameters). The program automatically resolves these dependencies and writes all internal parameters in terms of external parameters. The numerical values are required for the numerical evaluation of the analytically obtained results for λ hhh and are by default read from the UFO model. Moreover, one can change the default parameter values individually or all at once by specifying e.g. a SLHA [75,76] input file (see Section 4.3 for examples).
It should be stressed that the relations between internal and external parameters given via the UFO model in A.4) are essential for the renormalisation procedure. In particular Eq. (3) will be evaluated once all parameter dependencies have been applied. Thus, if any of the Higgs masses is not defined as an input parameter in the UFO model, it cannot be renormalised in the OS scheme automatically. Instead, the corresponding counterterm contribution would need to be provided manually via the custom counterterm δ (1) custom-CT λ hhh . It is, therefore, recommended to align the chosen parametrisation for input parameters in the UFO model along with the chosen renormalisation schemes.
The following steps are performed automatically using the information gathered before: • B.1) The UFO model is loaded and several checks are performed: -Whether all relevant couplings are present (especially quartic scalar couplings which are sometimes excluded in UFO outputs).
-Whether all relevant couplings are defined through the same Lorentz structures that are also used by anyBSM. Otherwise, one can use the model converter discussed in Appendix A.1.
• B.2) Definition of the SM-like particles and parameters based on the inputs made in A.2) in the file schemes.yml. The program is also capable of finding the SM particles and parameters automatically based on their PDG codes and numerical (mass) values. This functionality is also used to cross-check the user-input in order to avoid erroneous configurations.
• B.3) and B.4) All possible field-insertions into the generic diagrams are determined. The corresponding couplings and masses are inserted into the generic results. The calculation of any n-point function involved in the counterterm contributions follows the same procedure. Finally, the result for every n-point function is stored on-disk for caching/later use (see Appendix E).
• B.5) Collection of the individual results and construction of the expression for the renormalised λ hhh .
• B.6) Numerical or analytical evaluation. For diagrams with non-zero external momentum, the loop functions are evaluated using pyCollier (see Appendix F), which is a Python interface for COLLIER [77]. The analytical evaluation can be simplified/modified using sympy [78]. We want to stress that this particular strategy for obtaining a prediction for λ hhh in a given model has a number of ingredients in common with the calculation of many other observables. For this reason, the code anyH3 for calculating λ hhh is embedded into a larger program called anyBSM. The program anyBSM provides many utilities capable of performing the steps described above to set up the calculation of a particular observable. Utilities of this kind are for instance the interface to UFO or the insertion of UFO particles and Feynman rules into generic Feynman diagrams. These ingredients are used in submodules -of which anyH3 is the first one that has been implemented -to define actual quantities to compute. In fact, anyH3 only takes care of B.5) of the program points mentioned above while all other steps are taken care of in decoupled classes/modules of the program anyBSM. This class and module structure is depicted in Fig. 2. The class anyModel encodes all information about the used model (Feynman rules, particles, interactions, etc.). This information is then inherited by the anyProcess class which is used to calculate generic quantities like two-or three-point functions. This class makes use of various internal and external modules to derive the necessary diagrams and loop functions. The generic results of the anyProcess class are then used by the anyH3 class to specifically calculate the trilinear Higgs coupling. Future extensions of anyBSM featuring the calculation of new observables can be built on the basis of the anyProcess class (indicated by the ellipsis). In the end, the anyBSM class collects the classes for the different observables into a single object.

Program Tutorial
In this Section, we describe the basic features of anyBSM and anyH3. As explained above anyH3 is part of anyBSM, which provides a flexible Python framework for precision calculations (not only of λ hhh ). This Section is not meant to be a detailed manual but instead is intended to give an overview of the overall functionality. A detailed description of all available methods and options is available in the online manual, which can be found at https://anybsm.gitlab.io/anybsm.
Running the code requires at least Python version 3.5. The code is most easily used by installing the corresponding Python package by running pip install anyBSM which will automatically download and install anyBSM as well as all necessary dependencies. One necessary requirement not handled automatically by pip is the presence of a Fortran compiler (required for the compilation of COLLIER for the pyCollier dependency) such as gfortran or CLANG which can be installed from the systems package repository. Upon the first run of anyBSM, the model repository will be downloaded and saved into a user-specified location. 5 . The model repository is available online at https://gitlab.com/anybsm/anybsm_models and can be updated using the version control system git. The model repository will be expanded over time, and community contributions, in particular git merge requests, for tested models are welcome.

Basic syntax
anyBSM can either be integrated into Python scripts as a Python package or run directly from the command line. In addition, a Mathematica interface exists as well.

Python package mode
After starting a new Python session and importing anyBSM via from anyBSM import anyBSM a model -here for instance the SM -can be initialised via SM = anyBSM('SM') 5 The path to this location is defined in the config file which can be found at~/.config/anyBSM/anyBSM_config.yaml (Linux) or /Library/Preferences/anyBSM/anyBSM_config.yaml (Mac OS). By default, it is set to the folder models at the same location as the config file (e.g.,~/.config/anyBSM/models for Linux.
Alternatively, also a path to a UFO model directory can be given. As an overview, the dictionary anyBSM.built_in_models contains a list of all pre-installed UFO models and their installation directories. During the initialisation step, anyBSM will try to automatically identify the SM-like Higgs boson, for which λ hhh is calculated, as well as all other SM particles.
After the model initialisation, λ hhh can be calculated by running Here, "total" denotes the total value for λ hhh in GeV; "treelevel", the tree-level value; "genuine", the genuine one-loop contribution; "wfr", the contribution from external-leg corrections (diagonal and off-diagonal); "tads", the contribution from tadpole diagrams; "massren", the contribution from mass renormalisation; "vevren", the contribution from the renormalisation of the electroweak vev; and, "customren", the contribution from a custom counterterm. It should be noted that by default λ hhh is evaluated for vanishing external momenta. The option for using non-zero external momenta is described in Section 4.3.

Command line mode
As an alternative to using anyBSM as a Python package, it can also be called directly from the command line. It should be noted that the command line tool provides access only to the basic functionalities of the anyBSM library, unlike the Python package mode and the Mathematica mode discussed below.

Mathematica mode
The Mathematica interface can be conveniently installed as follows: Import["https://gitlab.com/anybsm/anybsm/-/raw/main/install.m"] InstallAnyBSM[] which checks for all requirements and adds the anyBSM interface to Mathematica's $Path variable. Afterwards, the interface can be used as follows: The result lambda is a Mathematica Association object similar to the Python dictionary object obtained in Section 4.2.1 using the Python library. However, by default the analytical rather than numerical results are returned, after conversion to valid Mathematica expressions. A list of all available functions (such as for e.g. the calculation of self-energies and tadpoles) within Mathematica is stored in the variable $AnyFunctions. More information is given in the online documentation. Furthermore, comprehensive Mathematica notebooks that demonstrate the use of anyBSM's Mathematica mode are provided in the examples repository.
The Mathematica mode has access to the full functionalities of the Python backend (i.e. the anyBSM library). For the sake of clarity, we will restrict ourselves to a description of the Python package in the next sections.

Setting parameters
anyBSM uses the default parameters defined in the respective UFO model. To change e.g. the value of the top-quark mass in the example SM calculation discussed above, we can run where in this example an effective (running) top-quark mass of 165 GeV is used. Alternatively, a LHA file can be used as input,

SM.setparameters('/path/to/LHA_file')
The program also defines a few additional UFO parameters in case they are not found in the UFO model. For instance, if the model does not define an external parameter named Qren (used for the renormalisation scale Q ren. ), the code introduces it internally with a default value of Qren = 172.5 GeV. The full list of additionally introduced parameters is discussed in Appendix D as well as in the online documentation.
The external momenta entering the computation of λ hhh can be specified by passing the momenta attribute to the lambdahhh function, e.g.

Renormalisation
Information about the renormalisation is saved in the file schemes.yml in the model directory. A simple example file for the SM is Here, the first SM_names block defines the names of various SM fields and parameters. The renormalization_schemes block can be used to define different renormalisation schemes. In the present example, the scheme OS is defined such that the mass of the field h as well as the VEV counterterm are renormalised in the OS scheme. This scheme is set as default scheme via the default_scheme directive. In addition, the scheme MS is defined so that the mass of the field h as well as the VEV counterterm are renormalised in the MS scheme. It should be stressed that this does not mean that all inputs are converted from OS to MS parameters but rather that the physical interpretation of these parameters is changed from OS to MS. However, for a consistent conversion of the parameters, all ingredients (i.e. twopoint functions) are provided by the program. A proper conversion between the schemes will be demonstrated in Section 6.1.
If a non-default scheme should be used, this can e.g. be specified during the model initialisation: As an alternative to using schemes predefined in schemes.yml, renormalisation schemes can also be generated interactively during the run time by using a new name that is not yet used in the schemes.yml file for the scheme directive during the model initialisation or by calling e.g. SM.add_renormalization_scheme('MS') afterwards. The new scheme will then be saved into the schemes.yml file. It is also possible to change the renormalisation scheme, e.g. between two calls of SM.lambdahhh(), using the appropriate method:

SM.load_renormalization_scheme('OS')
If no schemes.yml file is present in the UFO model directory, it will be generated automatically upon the first creation of a renormalisation scheme which automatically searches for all SMlike parameters (particles) based on their numerical (mass) values (and PDG identifiers), cf. Section 3.

Evaluation modes and output formats
anyBSM supports three different evaluation modes: • abbreviations: all results are given in analytical form using the UFO coupling abbreviations (GC_1, GC_2, etc.); • analytical: all results are given in analytical form using the full analytical form for all couplings; • numerical: the numerical values for all parameters are inserted, and a numerical result is returned.
The evaluation mode can be set e.g. via

SM.set_evaluation_mode('analytical')
if using anyBSM as a Python package. The default evaluation mode is numerical.
A detailed breakdown of the results (including results for individual diagrams) in the form of a PDF document can be produced by using draw = True as an additional argument for the lambdahhh function or via the -t option when using the command line interface

SM.lambdahhh(draw = True)
The individual results listed along with the diagrams are represented in a way which depends on the chosen evaluation mode (e.g. numerical or analytical/using abbreviations). The resulting PDF file is saved to the current working directory as well as the model directory. In order to make use of this feature, L A T E X needs to be installed.
In addition to the Mathematica mode, analytical expressions can be exported from within a Python session to Mathematica with the help of SymPy from sympy import mathematica_code mathematica_code(<sympy_expression>) Note that anyBSM includes a caching system which automatically saves the analytic results into json files (into the cache directory located in the model directory). This leads to a significant speed-up of consecutive runs, see Appendix E.

Getting Help
All Python classes and methods defined in anyBSM and anyH3 have meaningful doc-strings which can be issued by e.g.
from anyBSM import anyH3 help(anyH3) help(anyH3.lambdahhh) or directly using existing class instances (such as help(SM.lambdahhh) in the examples above). In addition, the online documentation makes use of these doc-strings and provides a search functionality.
The usage of the command line tool anyBSM is returned by the command anyBSM -h. For a given model, one can obtain further help by issuing the command anyBSM <model name> -h from the command line. The Mathematica interface of anyBSM also provides documentation for all its functions by issuing ?<function name> such as e.g. ?lambdahhh. Furthermore, it provides a list of available functions stored in the variable $AnyFunctions.
In addition, the anyBSM examples repository provides basic and concrete examples for all three interfaces and for the generation of new model files.

Built-in models and cross-checks
The models currently distributed alongside anyH3 are • the Next-to-Two-Higgs-Doublet Model (NTHDM)i.e. the real-singlet extension of the THDM; • triplet extensions of the SM with either a real triplet with hypercharge Y = 0 or a complex triplet with Y = 1. The two theories are denoted respectively TSM Y =0 and TSM Y =1 ; • the Georgi-Machacek model (a general version, as well as an aligned version); • a U (1) B−L extension of the SM (BmLSM); • the Minimal Supersymmetric Standard Model (MSSM).
These models, and associated conventions, are described in more detail in Appendix C. We emphasise again that additional models can also be included by the user in a convenient and fast way, as described in Appendix B. In this section, we present details about a variety of analytical and numerical cross-checks we performed to validate anyH3.

Cross-checks using analytical computations
To cross-check the routines implemented in anyH3, we compared the analytical results to calculations performed using FeynArts and FormCalc. We found full agreement for the following pieces of the calculation performed in anyH3: • Higgs and Goldstone boson self-energies (and momentum derivatives thereof) including self-energies with two distinct external scalars as well as charged scalar self-energies; • Higgs tadpoles; • vector boson self-energies including mixing self-energies (e.g. γ − Z mixing); • genuine one-loop corrections to scalar three-point functions; • one-particle-reducible contributions to scalar three-point functions.
These checks have been performed in the SM, meaning that all contributions to the renormalised trilinear Higgs coupling arising in the SM have been cross-checked. Contributions that do not exist in the SM (e.g. scalar-mixing self-energies) have been cross-checked in the THDM. As an additional cross-check, we have verified the cancellation of ultraviolet divergencies in the SM, the THDM, the SSM, the TSM Y =0 , and the TSM Y =1 . Moreover, we have found full agreement for the overall one-loop result for λ hhh with independent calculations in the SM and the THDM, performed with FeynArts and FormCalc, as well as with the results for the TSM Y =1 from Ref. [47] (see Appendix C for a more detailed description of the models).

Numerical cross-checks
In addition to analytical cross-checks, we have performed a series of numerical cross-checks by reproducing results from the literature.

SSM
As a first check, we reproduced the SSM results for κ λ derived in Ref. [38] (following the choice made in this reference and for the sake of comparison, we set the mixing between the CP-even states to zero). This reproduction is shown in Fig. 3 which is to be compared with Fig. 6 and 7 (upper-right) of Ref. [38]. In this figure, the momentum of two of the three external Higgs boson legs is always set on-shell p 2 1 = p 2 2 = m h 1 = 125 GeV. In the left plot of Fig. 3, the momentum of the third external Higgs leg is fixed at p 2 3 ≡ p 2 = 251 GeV and κ λ is shown as a function of the singlet mass. In the right plot, the singlet mass is fixed to 200 GeV and the external momentum of the third Higgs boson leg is varied.
The behaviour of both plots reproduces the behaviour found in Ref. [38]. However, the exact numerical values in the left panel are shifted due to different treatments of the externalleg corrections. The different treatments of external momenta also lead to a slightly different peak structure in the right plot. However, at p 2 = m 2 h 1 the different treatments coincide. To show this equality we use the shift to κ λ caused by the BSM sector, δ which was introduced in equation (25) of Ref. [38]. The external leg contribution to δ (1) hhh in the two different treatments reads Ref. [38] ext. leg is the one-loop Passarino-Veltman two-point function [79,80]. Thus, the two approaches yield the same external leg correction factors in the limit of p 2 → m 2 h 1 . We made use of this relation to cross-check the full analytical result obtained  with anyH3 with the result derived in Ref. [38] and found full agreement at p 2 = m 2 h 1 . For demonstrative purposes, we provide this cross-check using the Mathematica interface of anyBSM (cf. Section 4) as an example usage in the anyBSM examples repository.

TSM Y =1
As a further verification of anyH3, we reproduced results in the literature for the TSM Y =1 model [45]. Fig. 11 of Ref. [45] shows deviations of λ hhh from the SM prediction in the plane of the coupling λ 4 and the mass difference between the lightest and second-lightest BSM states (see Appendix C.7 for further details about the model). Our reproduction of this Figure is shown in Fig. 4. In the left panel, the lightest BSM states are the doubly-charged Higgs bosons; in the right panel, the lightest BSM states are the two neutral BSM Higgs bosons. Overall, we observe a very good agreement between our results and the results presented in Ref. [45]. The remaining small differences can be traced back to different SM input parameters used in Ref. [45].

MSSM
As a cross-check of the MSSM implementation, we reproduced the results of Ref. [31]. In this work, the leading O(m 4 t ) corrections to the trilinear Higgs coupling originating from scalar top quarks were calculated in the limit of vanishing electroweak gauge couplings. Setting the SUSY (and SM) parameters as in Ref. [31] (i.e., setting MQ = MŨ = 15 TeV, µ = |A t | = 1.5 TeV), we find very good agreement with their results (see our

Recovering the SM result in the decoupling limit
As an additional non-trivial cross-check of anyH3 (and also of the model files distributed alongside it), we have verified that the BSM contributions decouple if the masses of the BSM scalars are increased in a uniform way (see below for details), so that the SM result for λ hhh is recovered in this limit.
This verification is shown in Fig. 6, where κ λ is displayed as a function of the BSM mass scale M BSM . All BSM masses in each model have been chosen to be degenerate with each other with the mass value M BSM . The results for κ λ are shown in the SSM (green), the IDM (light blue), the THDM-II (dark blue), the NTHDM-II (stars), the TSM Y =0 (red), the TSM Y =1 (orange), and the Georgi-Machacek model (brown). In all models we chose appropriate input parameters such that the lowest-order couplings of the SM-like Higgs boson to the other SM states are exactly as in the SM. This setting is referred to as aligned scenario in the following. Further details about the chosen parameters of the models are specified in the caption. As a reference, the SM result for κ λ is indicated as a black line. It is clearly visible that for increasing M BSM κ λ quickly approaches the SM result for the chosen parameter settings. For M BSM 1 TeV, the deviations from the one-loop SM result are below ∼ 0.05. We note that, as will be seen in Section 6.4, values below the SM prediction are also possible in BSM models. Additionally, the details of the decoupling patterns in the different models strongly depend on the chosen parameters. Therefore, from the shown example no general conclusions can be drawn about how quickly decoupling occurs with increasing BSM mass scale in the various models.

Example applications
After having discussed cross-checks for validation, we present here a series of example applications. We first discuss estimates of the remaining theoretical uncertainties and then provide examples of results that go beyond existing studies in the literature.

Estimation of uncertainties in the computation of λ hhh
When comparing predictions for an observable (or a pseudo-observable) with experimental measurements or limits, an important consideration to ensure the reliability of the comparison is to estimate the theoretical uncertainty associated with the obtained prediction. We devote this section to discussing different contributions to the theoretical uncertainty associated with computations of the trilinear Higgs coupling: on the one hand, uncertainties due to missing higher-order terms, and on the other hand, parametric uncertainties due to the limited precision with which quantities entering the calculation of λ hhh are known experimentally.

Uncertainty from missing higher-order corrections in the SM
Focusing at first on the calculation of λ hhh in the SM, we begin by investigating the possible size of missing higher-order contributions. While in the SM two-loop corrections of O(α t α s ) and O(α 2 t ) are known [10,11,52], we refer here to all higher-order corrections that go beyond the full one-loop result that is obtained with anyH3 (see below for a discussion of the impact of the known two-loop corrections). Since calculations that are carried out at a given order in different renormalisation schemes differ by contributions that go beyond the calculated order, the comparison of different results of this kind for the same parameter point can be used as an estimate for the size of missing higher-order corrections -provided that the perturbative behaviour in the different schemes is of similar quality. We employ this method and compare the result obtained with anyH3 in the OS scheme with the one in the MS scheme.
In the SM, we use for the quantities that enter the tree-level expression of λ hhh -i.e., the Higgs-boson mass, the W -and Z-boson masses, and the electromagnetic coupling α em (the latter three quantities are in turn used to compute the Higgs VEV, see the discussion in Appendix C. 1

.2) -the following OS input values
where the notation M i indicates the OS mass of particle i (we deviate here from the lowercase notation m h employed for the Higgs mass in the rest of the paper in order to avoid ambiguities between OS and MS masses). This yields for the tree-level prediction of the trilinear coupling a value of λ In order to compare these values with those in the MS scheme, we must first convert the OS input parameters M h , M W , M Z , α em to the MS scheme. Working at Q = 172.5 GeV, and employing again an OS renormalisation of the tadpoles (see the discussion in Appendix C.1.1 for further details), we find after a one-loop conversion Using now these values as inputs for the calculation of λ hhh in the MS scheme, we obtain at Q = 172.5 GeV (again at full one-loop order, and with the same two choices of external momenta as above) The difference of about 0.3 − 0.4 GeV between the results obtained in the OS and the MS schemes constitutes a first estimate of a part of the unknown higher-order corrections; in relative size, the obtained shifts correspond to a difference of less than 0.2%. We note that if we had chosen to convert also the value of the squared Higgs mass used for the external momenta, the result for λ would have decreased by 0.2 GeV, giving rise to only a slight change of our uncertainty estimate.
Concerning the interpretation of the uncertainty estimates obtained so far, it should be emphasised that the scheme comparison done above in fact does not capture corrections to λ hhh involving the strong coupling α s , because the performed scheme conversions of the quantities M h , M W , M Z , α em do not involve this coupling at one-loop order. Such effects can on the other hand be estimated by converting the input value used for the top-quark mass (in contrast to the quantities entering the prediction of the trilinear Higgs coupling at tree-level, converting the top-quark mass entering at the one-loop level from the OS to the MS scheme directly gives rise to a two-loop effect in the prediction for λ hhh ). Starting from the OS value of M t = 172.5 GeV, a conversion including the leading O(α s ) and O(α t ) contributions to the top-quark self-energy (involving the strong gauge and top Yukawa couplings) yields an MS value of m MS t (Q = 172.5 GeV) = 166.3 GeV. Using this MS value in the computation of λ hhh , we obtain which corresponds to a total deviation of about 2 GeV compared to the OS results above where M t was used as input value. For the sake of comparison, we note that employing the OS computation of λ hhh using M h , M W , M Z , and α −1 em (0) from Eq. (7) but with the MS value of m t , we obtain As explained above, the dominant two-loop O(α t α s ) and O(α 2 t ) corrections to λ hhh are in fact known [10,11,52], and amount to about +3 GeV. Thus, our above estimates of higher-order corrections that are not included in the computation of anyH3 turn out to be close to the actual size of the known higher-order corrections. While for the case of the SM those two-loop corrections could be incorporated into the anyH3 prediction, in many other models for which anyH3 can be employed the corresponding corrections are not fully known. For reasons of uniformity we also restrict the SM prediction for λ hhh in anyH3 to the full one-loop level. An extension of the code providing the incorporation of higher-order contributions is left for future work.

Parametric uncertainties
Another source of theoretical uncertainty in the prediction of λ hhh arises from the experimental errors of the input parameters. In order to investigate the impact of these parametric uncertainties, we take into account the 1 σ ranges of the experimental input parameters as given in Ref. [81], while for the top-quark mass we use It should be noted that the variation of ±1 GeV for M t is indicated for illustration. The parametric uncertainty of the top-quark mass receives a contribution both from the experimental error of the measured mass parameter of ±0.30 GeV at the 1 σ level [81] and from the systematic uncertainty arising from relating the measured quantity to a theoretically well-defined top-quark mass. The parametric uncertainties that are induced by the masses of the other quarks and the leptons are negligible. We furthermore note that we do not consider a parametric uncertainty from the strong gauge coupling because it does not enter the expression of λ hhh at the one-loop level.
Varying each of the indicated experimental errors independently, we find the theoretical uncertainties induced in λ hhh shown in Table 1. As expected, the largest effect on λ hhh , with an induced uncertainty of ±0.5 GeV, originates from the experimental error of the mass of Parameter Exp. uncertainty Impact on λ hhh Table 1: Theoretical uncertainties in the calculation of λ hhh arising from the experimental errors of the input parameters.
the detected Higgs boson. Indeed the Higgs-boson mass enters the prediction for λ hhh already at the tree level, and while it is already known to a high level of accuracy its experimental error is still larger -by more than an order of magnitude -than the experimental errors of M W and M Z (and much larger than the parametric uncertainty associated with α em ). The theoretical uncertainties that are induced by the gauge-boson masses have only effects at the level of some tens of MeV or less. On the other hand, the experimental uncertainty of the top-quark mass has a stronger impact on λ hhh even though it only enters at the one-loop level. It should be noted that if in the future BSM parameters are measured, the parametric uncertainties in the prediction for for λ hhh induced by their experimental errors should also be taken into account.

Uncertainty from missing higher-order corrections in a BSM model: the example of the IDM
When considering the computation of λ hhh in BSM theories, BSM parameters can enter the tree-level expressions. This is not the case in aligned scenarios, like the IDM, where the treelevel prediction for the trilinear coupling is the same as in the SM and only the mass of the detected Higgs boson and the associated vacuum expectation value enter the lowest-order prediction. It should be noted that also in this case a comparison between the OS and MS results for λ hhh computed at one-loop order with anyH3 with a one-loop conversion of M h and v would not be sensitive to the type of contributions that can give rise to the largest effects at the two-loop level. For the specific case of the IDM one can infer from simple arguments of dimensional analysis that the leading two-loop corrections to λ hhh are of where Φ denotes either of the BSM scalars of the IDM, g hΦΦ is a coupling between the Higgs boson at 125 GeV and two BSM scalars, and λ 2 is the Lagrangian self-coupling of the inert doublet (c.f. Appendix C.4 for more details). These types of contributions are not generated by a one-loop conversion of M h or v.
Instead, the size of these contributions can be estimated via a one-loop conversion of the BSM scalar masses, which affect the size of the dominant one-loop corrections to λ hhh of O(g 3 hΦΦ /M 2 Φ ). In the following, we investigate for four different example benchmark pointslabelled BP1, BP2, BP3, and BP4 (defined in Table 2) -the potential size of these leading two-loop effects. We choose BP1 and BP2 with small splittings between the BSM scalar Inputs

MS masses
anyH3 results masses and the BSM mass parameter µ 2 , so that the couplings g hΦΦ (which are proportional to the difference M 2 Φ − µ 2 2 ) remain small, while we choose larger splittings for BP3 and BP4. Additionally, we set λ 2 to zero in BP1 and BP3, in order to investigate only terms of the form O(g 5 hΦΦ /M 4 Φ ), while for BP2 and BP4 we set λ 2 = 2 to also include effects of . We present the results obtained with the code anyH3 for the conversions of the scalar masses and for λ hhh in Table 2 (note that for the computation of (λ (1) hhh ) OS and for the scheme conversion of the BSM scalar masses, the tadpole contributions are renormalised on-shell).
As could be expected, we find that the OS and MS results are in very good agreement -differing only by 1.4% and 2.5% for the two choices of λ 2 -for the scenarios with small mass splittings (BP1 and BP2). For BP3 and BP4 featuring larger splittings, the discrepancy between the two results increases to about 5 − 6%. This confirms the known fact that the inclusion of two-loop corrections to λ hhh is increasingly important for parameter regions with larger splittings between the different BSM masses. The relative size of the O(λ 2 g 3 hΦΦ /M 2 Φ ) pieces compared to the O(g 5 hΦΦ /M 4 Φ ) ones also decreases for larger mass splittings, which simply follows from the lower power dependence on g hΦΦ ∝ (M 2 Φ − µ 2 2 ).

Comparison of renormalisation scheme choices for the TSM Y =0
Choosing a suitable renormalisation scheme is a crucial step for the calculation of κ λ . We illustrate this in Figs. 7 and 8 for the Y = 0 triplet extension of the SM (see Appendix C.6 for details about the model and its implementation). In Fig. 7, κ λ is shown for different renormalisation schemes as a function of the quartic interaction between the SM Higgs doublet and the BSM triplet (λ HT ) with fixed M H + = 1 TeV and λ T = 1.5. In the left plot, the tadpoles are treated in the FJ prescription, and therefore enter both the calculation of λ hhh and the parameter conversion explicitly, while  in the right plot an OS renormalisation is used for the tadpoles. If the mass of the SM-like Higgs boson and the EW VEV are renormalised in the on-shell scheme (red curve, identical in both plots), the dependence of κ λ on λ HT is very small. This is expected since the BSM masses are chosen at the TeV scale implying that all BSM corrections should be small as a consequence of decoupling. If the mass of the SM-like Higgs boson and the EW VEV are renormalised in the MS scheme (blue curves), the result for κ λ depends strongly on the choice of the renormalisation scale, as well as on the chosen treatment of the tadpoles. For all curves, m OS h = 125.1 GeV and v OS 250.7 GeV are used as input which are then converted in the first step to the MS scheme. Then, these MS quantities are used to calculate κ λ . For the considered scenario, the conversion can lead to very large shifts between the OS and MS quantities if the renormalisation scale is not chosen appropriately. If for example µ R = m t is chosen (solid blue curve), we encounter artificially large corrections to κ λ for large positive λ HT when employing FJ tadpoles. We note that this is due exclusively to the impact of the tadpoles on the MS parameters obtained from OS inputs, because the tadpole contributions in the calculation of λ hhh itself are the same independently of the employed scheme, as shown in Appendix C.1.1. Moreover, for this scale choice the MS mass of the SM-like Higgs boson quickly becomes tachyonic for negative λ HT in the case with FJ tadpoles, and for positive λ HT in the case with OS tadpoles. Similar issues appear for the choice of µ R = M H + (blue dotted curve) for which the MS mass of the SM-like Higgs boson becomes tachyonic for λ HT 1 (λ HT −1/2) when using FJ (OS) tadpoles. For µ R = (m t + M H + )/2 (blue dashed curve), however, the corrections are quite well-behaved, and a result close to the OS curve is obtained -although in the case of OS tadpoles, the Higgs mass once again becomes tachyonic for λ HT 3. Overall, the choice of OS tadpoles leads to more moderate effects in κ λ , and appears to be (if implemented) a preferable option for calculations in which scheme conversions of parameters are performed.
In the left plot of Fig. 8, we next present results for κ λ using OS inputs for the Higgs mass and EW VEV (and OS-renormalised tadpoles), but with the charged Higgs boson mass in the OS (red curve) and the MS scheme (orange curve). As in Fig. 7, we compare here results for the same parameter point, defined by λ T = 1.5 and an OS charged Higgs mass of M OS H + = 500 GeV (which for the MS curve is converted to the MS scheme). Because M H + only enters the prediction for λ hhh at the one-loop order, the difference between the red and orange curves is formally of two-loop order and, as was discussed for the IDM in Section 6.1.3, can serve as an estimate of the size of the unknown higher-order contributions to λ hhh . We can observe, as expected, that the results in the two schemes remain very close for small values of |λ HT 2|, and only differ for larger couplings. Finally, additional results for κ λ , calculated in the OS scheme, are shown the (λ HT , M H + ) parameter plane (for fixed λ T = 1.5) in the right panel of Fig. 8. 6 It is clearly visible that large corrections to κ λ are obtained for low M H + and large |λ HT |. The solid black contour lines indicate the parameter region that is excluded by the current LHC bounds on κ λ [16], 7 while the dashed black contour lines show the region that will be probed with the projected sensitivity based on the full HL-LHC dataset [19].

Comparison of BSM effects arising from mass splittings
In Section 5.2.4 we showed for various models that the BSM contributions to λ hhh vanish once all BSM states are decoupled simultaneously in an appropriate way. This behaviour is in accordance with the decoupling theorem [82] which states that, in the case of heavy new physics, all BSM contributions can be incorporated into coefficients C (n) i of higherdimensional SM operators, such that all BSM effects vanish for M BSM → ∞. A crucial requirement for this decoupling behaviour is that the C (n) i are small and do not increase with M BSM . Mass splittings between the BSM particles can modify this behaviour. This can happen for instance if some of the masses of the heavy BSM states φ BSM are mostly generated via the comparably small SM VEV, where φ 2 BSM schematically stands for the quadratic term of some BSM scalar, and Φ SM is the SM-like doublet. Thus, φ BSM receives mass contributions both from the quartic coupling λ XH and the mass parameter M L . For the case of M L ∼ v a large scalar mass M BSM v can be realised via λ XH 1. The quartic interaction of BSM scalars to SM-like Higgs bosons can lead to large contributions in this case, Another way of understanding the origin of such large contributions is related to the symmetry argument of the decoupling theorem i.e. that the decoupling of a heavy particle must not break any symmetries of the resulting effective theory (EFT). To demonstrate this we consider the states X i of some irreducible SU (2) L multiplet X with masses M X i =L ≡ M BSM and M X L ≡ M L M BSM . Taking the limit M BSM → ∞ leads to an EFT which is not SU (2) L -invariant anymore as X L cannot be incorporated into a smaller SU (2) L multiplet. As a consequence, portal couplings as in Eq. (16) between the SM-and the BSM-Higgs bosons become large for a large mass splitting. In fact, for all models implemented in anyH3 with additional states charged under SU (2) L we found the couplings c hhX i X j to behave as in Eq. (17), provided that appropriate parameterisations for the input masses are assumed. Note that the above discussion applies not only to the SU (2) L gauge symmetry but more generally to global and gauge symmetries of the BSM theories.
It should be noted that there can be regions in parameter space where the splitting of the mass parameters M 2 BSM − M 2 L is relatively large compared to the electroweak scale while the model can still be described perturbatively. From the THDM it is known that such large couplings can be constrained by the current experimental bounds on κ λ while being in agreement with all other experimental and theoretical constraints [13]. With the help of anyH3 one can easily go beyond the THDM and study the effect of couplings determined by Eq. (17) onto κ λ in other SM extensions. To demonstrate this, we use the example of the different SU (2) L extensions from Section 5.2.4 and fix one of the BSM scalar masses to M L = 400 GeV rather than having all masses degenerate at M BSM . Fig. 9 shows the resulting κ λ prediction as a function of M BSM for the IDM (light blue), the THDM (blue), the TSM Y =1 (orange) and the Georgi-Machacek model (brown). We want to emphasise again that all shown parameter points are chosen to be in the alignment limit, i.e. they have a tree-level prediction of κ We explicitly verified, using the anyPerturbativeUnitarity module of anyBSM, that all models fulfil the tree-level perturbative unitarity constraint in the high-energy limit for all shown values of κ λ .
It is also important to stress that this discussion is not restricted to the SU (2) L gauge symmetry of the SM but also applies to any other symmetry within or beyond the SM. 8

New results for λ hhh in the NTHDM
In this section we discuss an example of new results for λ hhh in the NTHDM. The scalar sector of this model contains three CP-even scalars, which can mix, and in turn three mixing angles are required to diagonalise the 3×3 CP-even scalar mass matrix (see model definitions and details in Appendix C.5). It is of interest in this context to investigate the potential impact of mixing on the prediction for the trilinear Higgs coupling.
As a brief illustration of phenomenological studies made possible with anyBSM, we present in Fig. 10 results for κ λ as a function of the second CP-even mixing angle α 2 . We choose a scenario of the NTHDM where h 2 is identified with the detected Higgs boson at 125 GeV, so that the alignment limit is reached for α 1 + α 3 → β − π/2 and α 2 → π/2. We consider in Fig. 10 scenarios of the NTHDM and THDM were the BSM scalars (h 2 , h 3 , A, H ± for the NTHDM; h 2 , A, H ± for the THDM) are mass-degenerate with a mass value of 300 GeV. We fix the BSM mass scales of both models (μ for the NTHDM and M for the THDM) to 100 GeV, resulting in a sizeable mass splitting giving rise to a significant contribution to κ λ . Additionally, we set tan β = 2 and α 1 + α 3 = β − π/2, while for the singlet VEV, v S , we adopt two values: v S = 300 GeV (blue curves) and v S = 3 TeV (red curves). Regarding our choice of renormalisation scheme, we employ here an OS renormalisation of all scalar masses and of the Higgs VEV, while the other parameters are renormalised in the MS scheme. Note that this scenario is devised as a simple setting in which to demonstrate calculations that are made possible by anyBSM. The range of parameters shown in Fig. 10 is not allowed in its entirety, however, we do not explicitly indicate the exclusion limits since we are not aiming at a thorough phenomenological analysis here. The solid curves in Fig. 10 show the full one-loop results for κ λ , while the dashed lines correspond to the tree-level results. For α 2 → π/2, we observe -as expected -that we recover the alignment limit, and for both possible values of v S the tree-level and one-loop predictions for κ λ converge to the results in the THDM, indicated by the black horizontal lines. In this limit, the additional singlet decouples entirely, and the dependence of κ λ on v S vanishes. A sizeable BSM contribution remains in this limit, yielding a value of κ λ ∼ 1.25, which arises from the corrections involving the THDM-like scalars (h 2 , A, and H ± ). On the other hand, away from the alignment limit, and for α 2 π/4, the relative importance of the loop corrections to λ hhh decreases significantly. It should be pointed out here, for completeness, that deviations of κ t from the SM are already constrained by experimental data to be below O(20%) (see for instance Ref. [14]). This implies that values of α 2 π/4 are already excluded in this scenario.
Furthermore, we can observe the interesting feature that the prediction for κ λ becomes negative as α 2 decreasesi.e. as one departs from the alignment limit. At this point, it is however important to remark that the sign of λ hhh is not a physical observable. A quantity that is of physical relevance is the relative sign between the trilinear Higgs coupling and other couplings of the Higgs boson, e.g. its coupling to top quarks. For this reason, we also present in Fig. 10 results (green line) for the coupling modifier of the top Yukawa interaction, which we denote κ t , at tree level. We find that for the entire range of α 2 , κ (0) t remains positive, so that a change in the relative sign between the trilinear Higgs coupling and the top Yukawa does occur -this can in principle lead to a significant increase in the Higgs pair production cross-section, as the destructive interference between the box and triangle diagrams occurring in the SM is avoided in this case. We leave a more thorough investigation of scenarios with negative trilinear Higgs couplings in the NTHDM for future work.

Momentum-dependent effects in λ hhh
By default, anyH3 evaluates the trilinear Higgs coupling setting the momentum of all external legs to zero. The code, however, also allows the evaluation of λ hhh for finite momenta (via the argument momenta of the lambdahhh function).
We demonstrate this for the THDM of type I (see Appendix C.3) and for the Y = 1 triplet extension of the SM (TSM Y =1 , see Appendix C.7) in Figs. 11 to 13. In all three scenarios, κ λ = 1 at the tree level. For the THDM-I (Figs. 11 and 12), we present results for two scenarios, with M = m h 2 = 400 GeV, m A = m H ± = 700 GeV for the first and M = m h 2 = 600 GeV, m A = m H ± = 1000 GeV for the second, and with tan β = 2 for both. Next, for the TSM Y =1 (Fig. 13), we set m D ±± = 400 GeV, m D ± = 500 GeV, and λ 4 = 4. All three scenarios are chosen such that significant BSM effects occur in the trilinear Higgs coupling, and additionally the second THDM-I scenario is devised specifically to obtain a value of κ λ larger than the current upper experimental bound of 6.3.
In the upper panels of Figs  renormalisation, see also Section 2.3. This means that the result shown in blue corresponds precisely to the quantity that enters the evaluation of the triangle diagram contributing to di-Higgs production. Comparing the solid blue and the orange-dashed curves, we observe that, for low to intermediate ranges of p 2 , the momentum effects are small in comparison to the overall size of the BSM effects, which shift κ λ to about 2.95, 7.9, and 3.3 respectively for the three figures. It is only for larger values of p 2 -namely p 2 600 − 700 GeV for the THDM-I scenarios and p 2 1.2 TeV for the TSM Y =1 -that the momentum effects become sizeable and cause a significant decrease in κ λ . Finite external momenta can also induce imaginary parts for λ hhh (for the calculation of κ λ , we take the real part of λ hhh ). These are shown in the lower panels of Figs. 11 to 13. Several particle thresholds are visible (corresponding to what can be seen in the upper panels of the corresponding figures): e.g., the di-Higgs threshold around p 2 ∼ 250 GeV, the ditop threshold around p 2 ∼ 350 GeV, and for instance for the TSM Y =1 (Fig. 13), also the D ±± D ∓∓ threshold ( p 2 ∼ 800 GeV) and the D ± D ∓ threshold ( p 2 ∼ 1000 GeV). Note that for the THDM-I (Figs. 11 and 12), there are no h 2 h 2 thresholds, because the λ h 1 h 2 h 2 coupling vanishes (due to the equality M = m h 2 ) and hence diagrams with internal CP-even scalars h 2 do not contribute to λ hhh . In Fig. 11, the AA/H + H − threshold is visible, as expected, at p 2 = 1.4 TeV.
The results shown in Figs. 11 to 13 can directly be applied to di-Higgs boson production by treating λ hhh as a momentum-dependent quantity entering the cross-section calculation. In this context, it is important to note that the di-Higgs invariant mass distribution typically  peaks around p 2 ∼ 400 GeV (see e.g. Ref. [84]) and then quickly falls off (by several orders of magnitude) as p 2 increases. This implies that the sizeable momentum dependence found for larger values of p 2 only has a small impact on the total di-Higgs boson production cross section. For the representative value of p 2 = 400 GeV, i.e. around the peak of the di-Higgs differential cross-section, the momentum-dependence contributes positive shifts of 3.4% for the first THDM-I scenario, 1.0% for the second THDM-I case, and 3.6% for the TSM Y =1 scenario. Furthermore, we observe that, in all three considered scenarios, the imaginary part of λ hhh remains minute until p 2 350 − 400 GeV (i.e. until the di-top threshold), and only reach sizeable values for p 2 well above 400 GeV, and thus far from the peak of the di-Higgs invariant mass distribution. Consequently, we find that evaluating λ hhh at zero external momenta is a good approximation of the full result for the scenarios considered here. We leave a more detailed investigation of the impact of non-zero momenta in di-Higgs production for future work. Moreover, we observe that for the points with the most sizeable BSM deviations in κ λ -such as the second THDM-I scenario -the relative magnitude of the momentum-dependent effects are the smallest, compared to the overall value of κ λ . Similarly, the relative size of the imaginary part of λ hhh , compared to its real part, appears to diminish for scenarios with larger κ λ .

Non standard couplings
While Section 2 discusses specifically the calculation of the trilinear coupling of the SM-like Higgs boson, the program anyH3 is in general able to calculate any trilinear self-coupling λ h i h i h i . To demonstrate this feature we again consider the real singlet extension of the SM, the SSM, and compute both the SM-like Higgs coupling λ hhh and the singlet coupling λ sss . For simplicity, we set the singlet-doublet mixing angle to zero, α = 0, which leads to the tree-level expressions , and (18a) where we fix the mass of the SM-like Higgs boson to m h = 125 GeV and allow the singlet mass m s to be either larger or smaller than m h (see Appendix C.2 for more details about the model). In this case the renormalisation of λ hhh at the one-loop order is identical to the SM. For the trilinear singlet coupling, we choose to renormalise κ S and κ SH in the MS and m s in the OS scheme. The renormalisation of the singlet VEV v S involves contributions from the singlet tadpole t s . anyH3 provides the necessary ingredients to compute the one-loop prediction for λ sss in the considered scenario.  200 GeV. This corresponds to the scenario λ ΦS = 1 that has been studied in Ref. [38]. 9 The corrections to λ hhh reach up to 20% for small values of m s in this scenario. For increasing m s the prediction for λ hhh approaches the one of the one-loop SM result. We find corrections to λ sss between −50 to +25% in the considered singlet mass range. For non-zero soft-Z 2breaking parameters the result for λ sss depends on the chosen renormalisation scale (while λ hhh does not depend on the renormalisation scale since the MS parameters κ S and κ SH do not enter the prediction for λ hhh at the tree-level). In order to demonstrate the dependence on the renormalisation scale we calculate λ sss for the two options of using a fixed scale µ R = m t and a dynamical scale µ R = (m t + m s )/2. We find that the difference in the λ sss prediction for these two scale choices is at least about a factor 3 smaller than the overall size of the corrections for most of the considered m s range. We have explicitly verified (besides the UV-finiteness of the result obtained with anyH3) that λ sss is independent of the renormalisation scale for κ SH = κ S = 0.

Conclusions
Obtaining information about the trilinear Higgs coupling λ hhh is crucial for determining the shape of the Higgs potential and for gaining a better understanding of the nature of the electroweak phase transition. BSM contributions to λ hhh can be large even for cases where all the couplings of the detected Higgs boson to gauge bosons and fermions are very close to the SM predictions. Thus, the comparison of the theoretical predictions for λ hhh in different models -taking into account contributions at the quantum level -with the available experimental constraints on λ hhh plays an important role for discriminating between the SM and extensions or alternatives of it.
It is therefore the main purpose of the public Python code anyH3, which we have presented in this paper, to provide precise predictions for λ hhh in a wide variety of models. anyH3, which is part of the broader anyBSM framework, calculates the trilinear Higgs coupling in the SM and renormalisable BSM extensions of it at the one-loop level. For the model input, the code supports the widely used UFO format. This allows the user to easily extend the library of models shipped alongside the code. Already 14 models are provided in this library.
The code implements generic one-loop corrections which are mapped to the respective UFO model. For renormalisation, semi-automatic routines allow the user to easily implement different renormalisation schemes. Besides calculating the trilinear Higgs coupling, anyBSM also supports the calculation of other quantities like scalar and vector boson self-energies.
We have validated the results of anyH3 by explicit analytical cross-checks, various consistency checks (such as cancellation of UV divergencies and decoupling of BSM contributions), and numerical cross-checks against known results in the literature. Besides comparing with known results, we have also presented new results for various models. In this context we have investigated different aspects like renormalisation scheme dependence, non-decoupling effects, momentum dependence, and negative trilinear Higgs couplings.
anyH3 can be used in the form of a Python module, called from the command line, or accessed via a Mathematica interface. All output quantities can either be evaluated analytically or numerically. To evaluate the required loop functions, anyH3 employs a link to the COLLIER library, which is available as the independent Python module pyCollier.
Besides detailed examples (including scripts to reproduce all plots in this paper), we also provide an extensive online documentation.
The code base of anyBSM is not restricted to the calculation of λ hhh , but can be easily extended to support the calculation of other observables like trilinear Higgs couplings with different external scalars or electroweak precision observables. We leave this for future work.

A Definitions and conventions for the generic results implemented in anyBSM
This appendix discusses the conventions employed in the generic expressions within anyBSM (and anyH3) and how to convert the conventions of an already existing UFO model to match  Table 3: UFO Lorentz structures used in anyBSM for vertices involving scalars S i , fermions F i , vectors V µ i , and ghosts U i . The four-momentum p µ j i = 'P(j,i)' in the second and third columns is carried by the field with the label i in the first column. In this table, γ µ denotes Gamma matrices, g µν the metric tensor, and P L,R the left-/right-handed projectors defined as P L,R = (1 ∓ γ 5 )/2. these. For generic vertices, which are used to obtain general results for Feynman diagrams, the program closely follows the conventions introduced in SARAH/SPheno [86], which are equivalent to the definitions used by the option InsertionLevel->{Generic} of FeynArts [65]. The UFO format, however, is much more general in the sense that every vertex V defined by a list of fields consists of two concurrent lists: a set of arbitrary Lorentz structures L i and a set of arbitrary couplings C i . The decomposition of V = L · C, however, is not unique and can differ between different specific UFO models. The explicit choice of Lorentz structures L i (which in turn fixes C for a given generic V ) for all renormalisable couplings used within anyBSM is shown in Table 3. If one wants to use anyBSM with a UFO model which does not obey the conventions of Table 3, the command line tool anyBSM_import (cf. next section) can be used to re-write all vertices of a given UFO model in terms of the Lorentz structures in Table 3.

A.1 Conversion between different conventions
As an illustrative example of how to handle different conventions, we consider the UFO implementation of the Inert-Doublet-Model (IDM) that is available at the FeynRules webpage as well as the UFO implementation which is included in anyBSM and was generated with the help of SARAH. The coupling between the neutral Goldstone boson and the top/anti-top quarks is defined in the FeynRules model as . This is obviously not compatible with the anyBSM conventions, in which all fermion couplings are written in terms of left/right-handed projectors, and anyBSM therefore would expect the vertex in the form For a variety of models, we checked explicitly that the prediction for λ hhh obtained with anyH3 is numerically identical when using the UFO model generated with SARAH (which can directly be used with anyH3) as well as with FeynRules after the use of the importer. 10 For instance, the U (1) B−L extended SM from the FeynRules model database can be converted in the following way: We have verified that the result obtained for λ hhh with this "converted" UFO model agrees with the one obtained with a version of the UFO model that was generated with the help of SARAH and built into anyBSM (cf. Appendix C.9), after adjusting all input parameters. 11

A.2 Available topologies and generic diagrams
Diagrammatically, the calculation of λ hhh (see Eq. (2)) can be expressed in the form  (19) where the solid lines are meant to be placeholders that are populated with all possible field-insertions (spin 0, 1 /2 and 1) from a given UFO model upon running anyBSM. The modelspecific results are obtained in anyBSM by inserting the Feynman rules into the result of the corresponding generic diagram. It should be noted that the tadpole contribution appears in both δ (1) tadpoles λ hhh and δ (1) tad. WFR λ hhh , and that the tadpoles as well as the external-leg corrections can optionally be turned off separately and instead be included in δ  CT λ hhh , have been discussed in Section 2.1 and will be discussed in more detail in Appendices C.1.1 and C.1.2. One crucial input for the renormalisation are scalar one-point functions, as well as all bosonic two-point functions, where S can be any scalar, and X, Y any scalar or vector boson. The internal solid lines are analogous to those in Eq. (19). Thus, anyBSM is at the moment able to calculate the scalar one-, bosonic two-point, and scalar three-point functions of any renormalisable QFT.

B Implementing new models
In this Appendix, we briefly describe the different possibilities and general strategies to create anyBSM-compatible UFO models. For detailed instructions and examples we refer to the online documentation and the examples repository, respectively.
Models in anyBSM are fully described using the UFO standard [57,58]. The code requires the following UFO files to be present in the model directory: • particles.py -specifying the fields present in the model, • parameters.py -specifying all input parameters (called "external" parameters) and "internal" parameters which are computed in terms of the input parameters, • lorentz.py -specifying all Lorentz structures appearing in the model, • couplings.py -specifying all couplings in terms of the model parameters, • vertices.py -specifying all vertices along with their corresponding Lorentz structures and couplings.
Additional Python files which are normally distributed for UFO models (e.g. function_library.py) are directly incorporated into anyBSM and are therefore not needed/ignored, cf. Appendix D.
UFO models for anyBSM can be created using e.g. SARAH [67][68][69][70][71] or FeynRules [57,[87][88][89]. However, the anyBSM package uses a specific convention for the Lorentz structures that are used in the derivation of generic results, cf. Appendix A. These conventions are identical to those of SARAH. For UFO models from other sources it is, however, not guaranteed that all vertices are expressed in terms of the Lorentz structures used in anyBSM. In such cases the tool anyBSM_import, which is discussed in Appendix A.1, can be used to test a given UFO model for compatibility with anyBSM and to perform a conversion, if required.
Nevertheless, there are currently some limitations on the models that can be used with anyBSM. The model • must not contain any non-renormalisable couplings (i.e., they are ignored), • needs to obey the conventions for Lorentz structures listed in Appendix A (otherwise the model-converter should be used), and • should define all SM-like particles, either via their PDG code as defined in Refs. [81,90,91] or by defining at least one renormalisation scheme in the schemes.yml file (see Section 4.4).
Additional current limitations are: • no external coloured states for the calculation of λ h i h i h i are supported (i.e., the scalar state under consideration must not be charged under SU (3) C ).
• The automatic OS renormalisation of the electroweak VEV (VEV_counterterm: OS in the schemes.yml) is only implemented for models with an electroweak ρ parameter that does not differ from one at lowest order, see the discussion in Appendix C.1.2.
• For the caching (cf. Appendix E) to work, the couplings defined in the file couplings.py need to be in ascending order (GC_1, GC_2, . . . , GC_23, GC_24). Otherwise only the caching of insertions should be used (set anyBSM.caching = 1).
New models can be added simply by including the corresponding UFO model files (run through the converter described above if necessary) in the models repository of anyBSM. The location of this directory is specified in the anyBSM configuration file which is written at the first start of anyBSM after the installation, cf. Section 4.1. Both locations, the config-path and the models-path, can be issued as follows import anyBSM from os import path print('anyBSM config file: ', path.join( anyBSM.config.appdirs.user_config_dir('anyBSM'), 'anyBSM_config.yaml')) print('anyBSM model directory: ', anyBSM.anyBSM.models_dir) To get an up-to-date list of known models contained in the model directory one can issue anyBSM -l within the command-line. For a list of built-in models, one can also consult the online documentation or visit the model repository which contains more information about the specific models and their implementation details. We also note that all underlying SARAH or FeynRules models, which have been used to generate the built-in UFO models, are publicly available in the examples repository.
Alternatively, new models that are unknown to anyBSM (i.e. stored outside of the models directory) can be used by providing the absolute/relative path to the UFO directory rather than just providing the model name (which is determined by the name of the model directory) at the initialisation step or as the first argument of the command-line tool.
For more information on how to implement a new model, we refer to the online documentation.

C Models currently provided in anyH3
In the following we review the models that are shipped alongside anyBSM and anyH3 in the UFO format. However, it should be stressed that the program is not restricted to this set of BSM models (and the SM) but works with all UFO models that fulfill the requirements described in Appendix B.
This appendix is organised as follows: In the first part addressing the SM, Appendix C.1, we in particular discuss the renormalisation of the trilinear Higgs coupling. In this context we describe the renormalisation of the tadpole and the VEV as well as the OS renormalisation of the mass of the Higgs boson. An OS mass renormalisation is also used as the default scheme in most of the BSM models discussed below (explicit examples on how to implement other schemes will also be given). In Appendices C.1 to C.10 we briefly describe the models in terms of their Lagrangian densities and chosen parametrisations as well as necessary conditions to achieve the alignment limit (i.e., the limit in which the the tree-level couplings of the SM-like Higgs boson are identical to the respective SM couplings). More detailed information about the individual models can be found in the online documentation and in the references therein. We emphasise once again that throughout this appendix the quantity λ hhh refers to the trilinear Higgs coupling of the detected Higgs boson with a mass of about 125 GeV.

C.1 The Standard Model (SM)
The SM Lagrangian used to generate the UFO model with SARAH is given by: The Yukawa matrices Y i are assumed to be diagonal for simplicity (i.e. the CKM matrix is unity). Their entries are traded for the measured lepton and quark masses used as inputs in the UFO model. We continue with the discussion of the tree-level scalar potential and its renormalisation. Assuming µ 2 < 0, the Higgs doublet of the SM obtains a VEV and is parametrised at the minimum of the potential in terms of the physical Higgs field h, the neutral (charged) Goldstone(s) G (G ± ) and the VEV v, Dropping all terms involving G ± or G, we can expand the tree-level potential as We choose to replace µ 2 and λ by defining the tree-level minimum t h as well as the squared tree-level mass m 2 h as follows In terms of m 2 h and t h , the potential becomes We introduce counterterms at one-loop order for the different parameters and the Higgs field entering the tree-level scalar potential as Turning now to the trilinear Higgs coupling, we find at the tree level that Correspondingly, the vertex counterterm (including field renormalisation) is given by At the one-loop order we use the parametrisation where δ (1) diag. λ hhh contains the one-loop diagrammatic corrections to the trilinear coupling (see the first two lines of Eq. (19)), and The second term in Eq. (32) arises from the counterterms that are associated with the diagrammatic tadpole contributions as described in Eq. (34) below. It turns out to have the same form as the tadpole counterterm contribution in Eq. (30) arising from the parametric dependence of λ where δ (1) tadpoles λ hhh refers to diagrams where a one-loop tadpole is attached to a Higgs-quartic coupling, i.e. δ (1) Rather than expressing λ hhh in terms of the tree-level Higgs-boson mass and the VEV, we want to express it in terms of physical inputs -namely the pole masses of the Higgs, W , and Z bosons (M h , M W and M Z ) as well as the fine-structure constant α em (0) (and ∆α, see Eq. (58) below). In practice, this can be done, either by performing conversions of mass parameters or by applying the OS scheme to fix the counterterms, with pole-mass relations of the form (using here the sign conventions of anyBSM) where we distinguish between a self-energy Σ no tad. , which does not contain one-loop tadpole insertions, and the full one-particle-irreducible (1PI) self-energy Σ. The transverse part of the gauge-boson self-energies is defined from the general decomposition Demanding that the tree-level input masses equal the pole masses, m 2 i = M 2 i , fixes the mass counterterms entering δ (1) CT λ hhh and δ (1) CT v in the OS scheme (see Appendix C.1.2 below). However, it is important to stress that, until this point, we have not yet specified the renormalisation of the tree-level minimum t h and of the electroweak VEV. In the following two sections we describe different possible treatments of the tadpoles and the VEV.

C.1.1 Equivalence of different tadpole renormalisation schemes
We now investigate different approaches for treating the tadpoles in the SM. The general tadpole contribution to λ hhh including all possible sources of tadpoles in Eq. (31) reads: The first term in Eq. (38) originates from the tree-level tadpole contribution, cf. Eq. (29), the second term arises from the vertex counterterm and the counterterms of the tadpole diagrams as described above, the third term contains the one-loop diagrammatic contributions to λ hhh of tadpoles (see the last diagram in the second line of Eq. (19)), cf. Eq. (34), and the fourth term arises from possible tadpole contributions in the mass counterterm. The nextto-last term denotes tadpole contributions to the VEV counterterm, while the last term vanishes since tadpole contributions to the field renormalisation (which is purely diagonal in the SM) drop out in the derivative w.r.t. the squared momentum. With Eq. (38) at hand, we can now discuss different choices of renormalisation schemes which amount to different relations/identities of/for the individual parts in Eq. (38). In the following, we restrict the discussion to the UV-finite parts of the various counterterms since UV-divergences are universal, and the UV-finiteness is not affected by the following discussions. 12 Tadpole-free MS scheme (TMS) A popular convention in e.g. Refs. [73,92,93] (this is for instance the default choice for loop calculations in SARAH/SPheno) is to employ MS renormalisation for the tadpoles, δ CT t h = 0, and to work at the minimum of the loopcorrected potential, which is realised by demanding that the total tadpole at one-loop order must vanish, i.e.
from which we can write where V (1) denotes the one-loop contributions to the effective potential. We note that this means that t h is formally of one-loop order. In addition, the VEV v is taken to be the true minimum of the loop-corrected potential. Using this in Eq. (38), we find at one-loop order the following total tadpole contribution to λ hhh , where we made use of Eq. (34) and Eq. (40) in the last step. The absence of the explicit dependence on any of the tadpoles is why this scheme is frequently called the tadpole-free MS renormalisation scheme (TMS). However, to compare the result in this prescription with the results obtained in a different renormalisation scheme, we need to express λ hhh (m 2 h ) -especially the tree-level piece -in terms of physical observables (i.e., λ hhh (M 2 h )), since m 2 h is not at the pole of the propagator in this scheme. The relation between the tree-level Higgs mass m 2 h and the pole mass M 2 h in this scheme reads Once inserted into the tree-level expression of λ hhh this gives an additional shift It should be noted that the counterterm contribution that is associated with the Higgs VEV, which is assumed to be the loop-corrected VEV in this scheme, does not explicitly introduce any tadpole contribution to λ hhh because no self-energy diagrams with tadpole insertions are included in this scheme.
OS tadpole renormalisation (tOS) An alternative treatment of the tadpoles was e.g. proposed in Refs. [94][95][96], where the one-loop tadpole corrections to the minimum of the tree-level potential (i.e., t h = 0) are required to be canceled by the tadpole counterterm such that the one-loop corrected minimum corresponds to the tree-level minimum (this is also often referred to as "OS tadpole condition" or "parameter renormalised tadpole scheme" (PRTS)), This has e.g. the advantage that no tadpole contributions explicitly contribute to the conversion between the MS and OS pole masses and that in general all diagrammatic tadpole contributions to any process are cancelled by the corresponding tadpole-counterterm diagrams (but tadpole counterterms furthermore appear in some other counterterms). Using these ingredients in Eq. (38) (together with Eq. (34)) we find which is identical to the result obtained in the TMS scheme, cf. Eq. (43). Note that we did not specify the scheme used for δ CT v. However, a cancellation of the one-loop genuine and counterterm tadpole contributions will occur separately therein due to the on-shell renormalisation of the tadpoles.
MS tadpole renormalisation at the tree-level minimum (FJ) Finally, we consider another scheme in which we again set t h = 0, but now renormalise the tadpoles in the MS scheme. Concretely, we require that δ CT t h cancels (only) the divergent part of δ (1) t h but is zero otherwise. This scheme is equivalent to the one known as the Fleischer-Jegerlehner (FJ) scheme [74]. Returning to the master expression, Eq. (38), of λ hhh we have in this scheme In this scheme we need to properly extract m 2 h as well as the VEV from their relation to physical observables including tadpole contributions (since they are not cancelled by their OS counterterms), which are then finally inserted into the tree-level expression for λ hhh . For simplicity, let us assume now that we extract the VEV from its relation to the OS pole mass of the Z boson M Z using MS values of the EW gauge couplings, 13 v(M Z ) ≡ 2 With the FJ treatment of tadpoles, this results in a tadpole contribution to the VEV counterterm that reads with This yields for the tadpole contribution from the VEV counterterm to the trilinear Higgs coupling In addition, we again need to make sure to express m 2 h in terms of the pole mass M 2 h , which we achieve by renormalising it on-shell. With MS-renormalised tadpoles, the finite shift relating the lowest-order Higgs boson mass to the pole mass contains a tadpole contribution of the form +3/vδ (1) t h . Thus we find a shift to λ hhh of the form +9/v 2 δ (1) t h .
Summing all contributions involving the tadpoles, we find a total of which is again in agreement with the two previous schemes.
In conclusion, the three tadpole schemes discussed here yield the same finite contribution (up to higher-orders) to the renormalised λ hhh of the form +3/v 2 δ (1) t h UV-finite , once we take into account that for a proper comparison we need to express all parameters in terms of physical observables. This result is expected since the relation between the input parameters that are expressed in terms of physical observables and the process of Higgs pair production corresponds to a relation between physical observables for which the actual treatment of the tadpoles must not matter.
Comparison of tadpole schemes and implementation in anyH3 Given the equivalence of the tadpole schemes demonstrated above, we are free to choose the most convenient treatment. Since the UFO format in general does not contain any information about the (tree-level) tadpoles, we choose the FJ treatment per default, since here we have t h = 0. While in the equivalence proof above we re-wrote all tadpole-inserted contributions in terms of tree-level couplings and propagators multiplied by the one-loop one-point function, in the actual code implementation we cannot automatically perform this re-organisation of the calculation but instead we generate and calculate all tadpole-inserted diagrams separately.
However, the aim of anyH3 is to be very flexible regarding the choice of renormalisation schemes. Therefore, the FJ scheme is only the default choice but the program easily allows the user to restrict to genuine loop diagrams and to add the tadpole contributions using a custom user-defined counterterm.
For example, in the SM the default FJ and the alternative tOS scheme can be defined in the following way in the schemes.yml file: where the "OS" scheme corresponds to the standard FJ scheme while the "OStadpoles" scheme corresponds to the tOS scheme and makes use of a custom counterterm "custom_CT_hhh" which is precisely the tadpole contribution that was derived above. 14 In addition "tadpoles: False" is used in the "OStadpoles" scheme to switch off all tadpole contributions in the calculation of the other counterterms, of self-energies, and of the loop corrections to λ hhh . The switch "VEV_counterterm: OS" refers to the calculation of the VEV counterterm contribution as described in the next section. It should be stressed that both, the "VEV_counterterm" and the "mass_counterterms" options, are compatible with the "tadpoles: False/True" option. In the example above this means that all 1PI self-energy diagrams, including tree-level propagators with one-loop tadpole insertions, are automatically taken into account in the scheme "OS" but not in the scheme "OStadpoles" where, because of the setting "tadpoles: False", Σ no tad. is used for all self-energies. Therefore, all tadpole contributions have to be added manually when using "tadpoles: False" via the "custom_CT_hhh" directive. The two different schemes can easily be compared numerically as follows: which shows perfect agreement within the numerical accuracy. We note that all other parameters in this example are renormalised in the OS scheme. If this is not the case, then a conversion has to be performed, resulting in numerical differences -see the discussion in Section 6.1. For numerical studies, a possible drawback of the FJ scheme is that is can suffer from numerical instabilities due to large numerical cancellations. In Appendix C.3, we show with a THDM example that the code is not affected by this issue in the phenomenologically relevant parameter region. However, as discussed in Section 6.2, if conversions of parameters are performed, we recommend to use an OS scheme for the tadpoles, whenever this is technically feasible, in view of the better behaviour of the perturbative series that this choice exhibits.

C.1.2 Vacuum expectation value renormalisation from OS quantities
In the previous discussion we used MS gauge couplings in the context of the extraction of the VEV within the FJ scheme. However, for practical reasons a direct determination of v in terms of measured quantities is desirable. Solving the tree-level relation between the electromagnetic charge e, the vector boson masses M Z/W , the weak mixing angle s θw ≡ sin θ w and the VEV v, for the electroweak VEV v, we can define the VEV in terms of the measured OS values of α em (0) (with e = 4πα em (0)), M W , and M Z as The counterterms for the parameters entering this relation are defined as where the term sng(s θw ) ensures that the sign convention in the covariant derivative is taken into account correctly -indeed a change in the sign with which the covariant derivatives are defined results in a sign flip in the sign of the weak mixing angle, and both sign conventions exist in the literature and can be employed when creating UFO model files. The transverse part of massive vector self-energies is defined in Eq. (37), while for the photon self-energy we have The resulting VEV counterterm in the OS scheme reads from which we can also relate v OS and v MS . In Eq. (54b), the photon vacuum polarisation Π γ (0) is split into three contributions: where Π γ (0)| heavy contains contributions from heavy fermions as well as all bosonic contributions and where ReΣ T γ (M 2 Z ) light is the transverse part of the photon self-energy considering only the light degrees of freedom of the SM, i.e. all leptons and the five light quarks. The contribution of the light fermions to the vacuum polarisation, Π γ (0)| light , would develop infra-red (IR) divergences in the limit of vanishing fermion masses. Thus, its contribution is absorbed in the quantity ∆α, defined as ∆α = ∆α (5) had. + ∆α lep. = 0.02766 + 0.031497687 , where ∆α (5) had is extracted experimentally [97], while ∆α lep. was computed in Ref. [98]. Note that this value can also be changed using the UFO format, cf. Appendix D. The evaluation of ReΣ T γ (M 2 Z ) light as well as of the heavy field contributions to the photon vacuum polarisation, Π γ (0)| heavy , is straightforward and numerically stable.
From Eq. (56), it can also be seen that the tadpole contribution in this treatment of the VEV is identical to the one chosen in the discussion of the FJ scheme in Appendix C.1.1: the tree-level couplings between the Higgs boson and the massive vector bosons normalised to their squared masses are universal (i.e. the same for W and Z) such that the tadpole contribution to δ (1) M 2 W /M 2 W and δ (1) M 2 Z /M 2 Z are identical. Thus the tadpole shift is only caused by vΣ T W W (M 2 W )/(2M 2 W ) in this scheme, which is identical to the one in Eq. (48) up to higher-order corrections.

C.2 The SM with a real singlet (SSM)
The most general potential that couples a real scalar gauge singlet S to the SM reads This model is implemented in the SARAH package under the name SSM. This implementation was used to create the UFO files after adapting the conventions for the input parameters of the SM sector, as well as those discussed in the following, to the anyH3 conventions.
After spontaneous symmetry breaking, the CP-even component of the Higgs doublet and the singlet are assumed to acquire VEVs, The fields s and h mix to two CP-even eigenstates h 1,2 with masses m h 1 < m h 2 : where c x ≡ cos x and s x ≡ sin x. We eliminate the parameters µ 2 and m 2 S using the tree-level tadpole equations which yields the following squared mass matrix: Furthermore, we choose to eliminate the three dimensionless parameters λ H , λ S , λ SH in favour of the two masses m 2 h 1,2 and the mixing angle α: Consequently, the scalar sector of the model is determined by the following input parameters: At tree level, the expressions of the trilinear self-couplings of the two CP-even states read In the limit κ S , κ SH → 0, the model obeys a spontaneously broken Z 2 symmetry and only adds three BSM parameters (the Z 2 breaking VEV, one mass as well as one mixing angle) to the SM parameters. For more detailed information about the model file generation we refer to the examples repository.

C.2.1 Alignment
The alignment limit, in which all tree-level couplings of h 1 to the SM sector (including the trilinear Higgs coupling) become identical to those of the SM, is achieved by choosing: However, in this limit the coupling λ SH is still non-zero and thus the one-loop prediction for λ hhh = λ (1) h i h i h i can differ from the SM prediction. Note that the SM-like Higgs boson can in principle be either h 1 or h 2 , see for instance the discussion of Fig. 15. In the following discussion we assume without loss of generality that h 1 = h is the SM-like and h 2 = s is the singlet-like Higgs boson.

C.2.2 Renormalisation
In the following we briefly describe the default renormalisation conditions in the (BSM) Higgs sector that have been implemented in the schemes.yml file. It should be emphasised that this choice is not fixed and can be changed by the user. Note also that in the following (as well as for the discussion of renormalisation in other models), we omit the superscript (1) and subscript CT on counterterms, as there should be no risk of confusion.  δv S : From the one-loop RGEs one can infer that the singlet VEV does not receive a UV-divergent counterterm at the one-loop order. As a simple cross-check of anyH3, we diagrammatically verified this fact. Therefore the singlet VEV counterterm does not introduce a renormalisation scale dependence in λ (1) h 1 h 1 h 1 . Thus, we could choose not to renormalise v S . However, due to the chosen default FJ renormalisation scheme in anyH3 for all appearing tadpoles (meaning that v S is the VEV of the tree-level potential) we can also consider a finite shift to the singlet VEV which ensures that it corresponds to the minimum of the loop-corrected potential: δα and δZ h 1 h 2 : The default choice in anyH3 is to calculate all external leg corrections at the given external momenta, which are set to zero by default. Thus, per default a fully OS wave-function renormalisation which properly removes the mixing between external states is not used. This is the default "OS" scheme (denoted "OS" because the SM-sector is still renormalised on-shell) defined in the schemes.yml file: In the example above the counterterm dvS corresponds to Eq. (68), while Derivative(lambdahhh_tree, 'vS') takes the derivative of the tree-level prediction for λ h 1 h 1 h 1 w.r.t. v S . The final result of the counterterm contribution again needs to be saved in the variable self.custom_CT_hhh, since this variable is internally used by anyH3. In this default scheme, α is renormalised in the MS scheme.
For demonstration purposes, we also implemented a non-minimal renormalisation scheme based on the OS scheme used in Ref. [99] for heavy Higgs decays. This version of the OS scheme fixes the mixing angle counterterm to The corresponding renormalisation scheme in the schemes.yml file is called "OSmixing": Note the additional entry wfrs: False which turns off the automatic calculation of all external-leg corrections (which defaults to p 2 = 0). Instead they have been added in the custom_CT_hhh-section using on-shell momenta for all three external legs. δκ S and δκ SH : The soft-Z 2 -breaking parameters are renormalised in the MS scheme. Thus, the trilinear Higgs self coupling of the SM-like Higgs boson is given purely in terms of OS parameters if we work in the Z 2 -symmetric limit where κ S , κ SH → 0.
We have explicitly verified that in this limit the result obtained in the "OSmixing" scheme is manifestly renormalisation-scale independent. We also provide a Mathematica notebook in the examples repository which explicitly demonstrates the UV finiteness of this scheme using the Mathematica interface of the anyBSM code.

C.3 The Two-Higgs-Doublet Model (THDM)
For the description of our implementation of the THDM we proceed along similar lines as for the SSM. The tree-level scalar potential of the built-in CP-conserving THDMs reads This potential obeys a Z 2 symmetry (Φ 1 → Φ 1 , Φ 2 → −Φ 2 ), which is softly broken by the m 2 12 term. Even though the different types in the Yukawa sector [100][101][102] play only a minor role in the prediction of λ hhh , we implemented for convenience four distinct models corresponding to type-I, type-II, type-X and type-Y (the latter two types are also often called type III and type IV, or flipped and lepton-specific) -see Ref. [103] for a review. After spontaneous symmetry breaking, the two doublets obtain the VEVs v 1,2 with v 2 /v 1 = tan β. The mixing between the two CP-even states is described by the parameter α. We eliminate µ 2 1 and µ 2 2 using the tadpole equations and we trade λ 1,2,3,4,5 for the masses of the two CP-even Higgs masses m h 1,2 , the CP-odd Higgs mass m A , the charged Higgs mass m H ± , as well as the mixing angle α. Thus, the chosen (BSM) input parameters are: For more detailed information we refer to the online documentation and the examples repository.

C.3.1 Alignment
The alignment limit of the THDM is achieved if the electroweak vacuum expectation value is aligned with the SM-like Higgs boson. If the lightest CP-even state is SM-like, this is achieved by setting sin(β − α) = 1; if the heavier CP-even state is SM-like, cos(β − α) = 1 corresponds to the alignment limit.

C.3.2 Renormalisation
The renormalisation schemes provided for the THDM are similar to the SSM. The default scheme only renormalises the SM sector OS, while for the mixing angles the MS scheme is used. As an additional option, which we employed for checking UV finiteness, we provide a scheme which uses an OS mixing angle as in the SSM, cf. Eq. (69). Furthermore, an OS counterterm for tan β can be either defined via the charged Higgs/Goldstone self-energies or via the self-energies of the pseudoscalar Higgs/Goldstone bosons. Since the two schemes differ by finite contributions, we have implemented both of them. The parameter M is always renormalised in the MS scheme.
A complete list of all renormalisation schemes defined in the schemes.yml file can be accessed by e.g. calling the list_renormalization_schemes() function: from anyBSM import anyBSM THDM = anyBSM('THDMII') THDM.list_renormalization_schemes() or by consulting the online documentation.
In the alignment limit (see Appendix C.3.1), the tree-level prediction for λ (0) hhh is identical to the SM (cf. Appendix C.1), which means that λ hhh can be fully expressed in terms of SM OS quantities in this case. Thus the tOS and the FJ schemes will formally yield the same result, as discussed in Appendix C.1.1. However, the FJ scheme is known to be perturbatively unstable due to large cancellations between different tadpole-inserted diagrams. To assess potential instabilities, we also implemented the tOS scheme in the alignment limit.
In the left panel of Fig. 16, we compare the numerical prediction for λ hhh in the THDM-II using the tOS and FJ schemes for a specific scenario with sin(β − α) = 1, tan β = 2, and m H ± = m A = m h 2 ≡ m Φ = √ M 2 + ∆ 2 (with ∆ = 200 GeV) as a function of M . In this parametrisation, it is expected that the BSM contributions decouple in the limit M v, where the BSM states become heavy and mass-degenerate, and that the THDM result converges towards the SM result (black line). This is indeed observed for both schemes near the scale M = 10 3 − 10 4 GeV. One can see that as long as M 10 5 GeV -well within the decoupling region -the tOS scheme (orange-dotted) agrees very well with the FJ scheme (blue-solid line). For M 10 5 GeV, numerical instabilities start to develop in the FJ scheme.
In the right panel of Fig. 16, we show the relative difference in λ hhh between the two schemes for various values of ∆. Even for very large couplings close to the unitarity limit (where the latter is reached for ∆ 400 GeV), implying very large loop corrections for

C.4 The Inert-Doublet Model (IDM)
The IDM is a variant of the THDM in which the Z 2 symmetry is imposed exactly. Its potential is for this reason identical to that of the THDM, with the exception that m 2 12 vanishes: After spontaneous symmetry breaking we only allow the SM-like Higgs doublet to obtain a VEV, Φ 1 = v / √ 2, but not the second doublet, Φ 2 = 0, because the Z 2 symmetry is exact. We again eliminate µ 2 1 using the tadpole equations as well as write λ 1,3,4,5 in terms of the scalar masses and µ 2 2 . Thus, the input parameters of the scalar potential are: with For more information about the model we refer to the online documentation.

C.4.1 Alignment
Since the CP-even Higgs bosons do not mix in this model, the tree-level couplings of the SM-like Higgs boson are automatically equal to the respective SM couplings.

C.4.2 Renormalisation
The tree-level prediction of λ hhh in this model is always identical to the SM. Hence, we only need to fix the renormalisation conditions for the SM parameters at the one-loop order as discussed in Appendix C.1.

C.5 The THDM with a real singlet (NTHDM)
The Higgs sector of the NTHDM consists out of two SU (2) L doublets Φ 1,2 and a real singlet S. In the model file, the discrete symmetries are imposed. Then, the most general potential reads 15 After spontaneous symmetry breaking, the CP-even components of the Higgs doublet ρ 1,2 and the singlet ρ S mix with each other. The mixing matrix is defined as where α 1,2,3 are the three mixing angles. Then,  The mixing angle of the CP-odd and charged components of the Higgs doublets is β (with tan β = v 2 /v 1 ), like in the THDM. The Yukawa sector of the model provided with anyH3 is defined as a type-II Yukawa sector. For more details about the NTHDM, see e.g. Ref. [104]. 15 We note that the off-diagonal mass term m 2 12 is defined here with the opposite sign compared to its counterpart in the THDM (the latter follows the convention in the SARAH model files).
We again demand that the triplet does not take part in electroweak symmetry breaking. Thus, there is no mixing in the CP-even sector. The mass parameter m 2 is eliminated using the SM-like tadpole condition. Then, the tree-level Higgs boson masses are given by We use m D ± , m D ±± , and the couplings λ 2,3,4 as BSM input parameters. All couplings of the SM-like Higgs boson in this model are as in the SM at tree-level. Therefore the renormalisation is treated as in the TSM Y =0 and the IDM. For more details about the Y = 1 inert triplet model, see e.g. Ref. [47].

C.8 The Georgi-Machacek Model (GM)
The Georgi-Machacek model is a triplet extension of the SM, which -in contrast to the previously discussed triplet extensions -respects the custodial symmetry of the SM at tree-level and therefore has a tree-level ρ parameter of one. The model adds one complex and one real triplet, respectively denoted ξ and η, to the SM scalar sector. The doublet and the triplets are effectively re-written in terms of one bi-doublet Φ and one bi-triplet X transforming under the custodial group SU (2) L × SU (2) R : The most general potential that is manifestly invariant under SU (2) L × SU (2) R reads [106] V GM = µ 2 2 which transforms the bi-triplet X in terms of Cartesian field-coordinates. The σ i (t i ) are the SU (2) generators in the adjoint (fundamental) representation. We provide two UFO models which implement the two scenarios where either M 2 X < 0 or M 2 X > 0.
In the first case the bi-triplet obtains a VEV v X that triggers mixing between the doublet and the triplet states, i.e., the SM-like Higgs will have a non-zero triplet component. We solve the tadpole equation for the SM-like doublet to eliminate µ 2 . The tadpole equation for the bi-triplet is used to eliminate M 1 rather than M 2 X , as which conveniently allows to take the alignment limit M 1 , v X → 0 (while keeping v X /M 1 finite). After spontaneous symmetry breaking we trade λ 1 , λ 2 and λ 5 for the SM-like Higgs mass, the triplet mass M 3 and the fiveplet mass as BSM input parameters. We do not trade λ 2 for m h 2 , but keep λ 2 as input since the parametrisation which uses m h 2 as input suffers from numerical instabilities in the decoupling limit (see also Ref. [107] for related discussions). For more information on the model we refer to e.g. Refs. [106,108] and the online documentation.

C.8.1 Alignment
Alignment in this model can only be achived with a vanishing triplet VEV (i.e., sin θ H = 0). Since this also significantly simplifies many tree-level relations, we provide a separate model (GeorgiMachacekAligned) without a triplet VEV (i.e., M 2 X > 0) in addition to the UFO model implementing the general (i.e. M 2 X < 0) case (GeorgiMachacek).

C.8.2 Renormalisation
For the GeorgiMachacekAligned UFO model the usual SM-like OS/MS renormalisation procedure is sufficient. In the general GeorgiMachacek-implementation all BSM parameters are assumed to be MS parameters for simplicity. Non-minimal renormalisation conditions in this model have been proposed e.g. in Ref. [46].
C.9 A simple U (1) B−L extension of the SM (BmLSM) We also provide a UFO model which extends the SM by a gauged U (1) B−L symmetry that is spontaneously broken by the VEV v X of a complex scalar S. The scalar potential reads We use the tadpole equations to eliminate µ and µ S as well as trade v x , λ 1 , λ 2 and λ 3 in favor of the Z mass M Z , the two scalar masses m h 1,2 , and the mixing angle α. Right-handed neutrinos are added in order to cancel the additional gauge anomalies: For simplicity we assume that the Yukawa matrices Y ν and Y x are diagonal. Thus, we can diagonalise the neutrino mass matrices, resulting in 6 eigenvalues which scale for v SM v X to first non-vanishing order as We choose to trade the components of Y v and Y x for the neutrino masses. In consequence, the BSM input parameters of the model are Kinetic mixing between the U (1) Y and U (1) B−L is not considered (i.e., set to zero at tree level) in this model implementation. For more information on the model see e.g. Ref. [109].

C.9.1 Alignment
All couplings of the SM-like Higgs boson to the SM sector are rescaled by cos α. Thus, α = 0 is the proper alignment limit in the Higgs sector. However, it should be noted that in this limit the bosonic BSM sector completely decouples from the SM and that the one-loop corrections to λ hhh are identical to the SM. Furthermore, in this limit charged lepton-currents are still not SM-like as they get modified by factors of O(m ν i /m ν i+3 ), i = 1, 2, 3. However, those contributions are numerically negligible.

C.9.2 Renormalisation
All SM parameters are renormalised using the automatic OS/MS procedures described above. The BSM sector is renormalised in the MS scheme. Alternatively, we provide an OS scheme for the Z mass as well as a scheme with an OS counterterm for α along the same lines as in the SSM and THDM.

C.10 Minimal Supersymmetry (MSSM)
A minimal supersymmetric version of the SM is also provided in terms of a UFO model. The UFO files were generated for the CP-conserving but non-minimal flavor-violating MSSM using SARAH. Therefore it is recommended to make use of the corresponding SPheno spectrum generator in order to calculate parameter points for this model (see e.g. the example in Section 6.6). An example SLHA file produced with SPheno is contained in the MSSM model directory and can be used e.g. in the following way:

C.10.1 Renormalisation
In contrast to the previously discussed models, the tree-level expression for λ (0) hhh in the MSSM does not depend on the Higgs boson mass, but on m Z , v SM , α, and tan β. We are again able to automatically renormalise all SM-parameters OS. tan β and α are renormalised in the DR scheme. Since anyBSM by default uses dimensional regularisation (i.e., the MS scheme), one needs to switch to dimensional reduction (i.e., the DR scheme) in the schemes.yml of the MSSM model file in the following way dimensional_reduction: True D Additions to the UFO standard While the UFO standard itself provides a very general strategy to store detailed information about a particular QFT, we found that several additions and adjustments in its actual implementation are useful in order to overcome difficulties in the renormalisation procedure as well as to improve computing performance.
As a first step, anyBSM expects various input parameter definitions to be present in the file parameters.py that are required for the renormalisation procedure. If the parameters do not exist in the UFO files, the program creates them using the UFO object library with the according default values: • Qren is a parameter used for setting the renormalisation scale and defaults to Q ren = 172.5 GeV (i.e. the top-quark mass).
• GFermi is used to set the Fermi constant. Alternative parameter names that are also searched for are Gfermi, gFermi, gF, GF and Gf. The default value is G F = 1.16637 · 10 −5 GeV −2 .
• Zsignfac is a parameter which fixes sng(sin θ w ) (i.e., the sign of the weak mixing angle, cf. Eq. (54b)). If it is not provided by the user, it is obtained automatically upon run-time using the method getSignSinThetaW() which determines the sign from the relative sign of the Z-top-top and photon-top-top couplings.
In addition, models are loaded with an anyBSM-custom object_library.py that adds some convenience and performance features: • The global all_<ufo object> variables are dictionaries instead of lists.
• Warnings appear if UFO objects with the same name (based on the .name-attribute) are initialised.
• A nvalue-attribute was introduced, containing the numerical value obtained with the current set of inputs.
• Similar to nvalue, a nmass-attribute for Particle-instead of Parameter/Couplingobjects was introduced. Note that a call of setparameters() (see Section 4.3) automatically updates the nmass and nvalue of all UFO objects.
• A new method UFOBaseClass.dump(), which returns a string representation to be interpreted by Python, was implemented.
• Particle.anti() now avoids creating duplicate Particle instances on each call.
• function_library.py was completely rewritten to avoid the call of __exec__ of functions. This yields a significant performance boost for large models of up to 3 orders of magnitude in run-time.
It should be noted that if any of the auxiliary files mentioned above is present in the UFO model directory, it is going to be ignored upon the model initialisation since only the modified anyBSM methods are used.

E Caching in anyBSM
One advantage of the anyBSM framework compared to many other similar tools is that it provides a sophisticated caching mechanism. This allows us to significantly improve the run-time of a phenomenological study. There are two levels of caching that can be set via e.g.
from anyBSM import anyBSM SM = anyBSM('SM', caching = <cache-level>) with <cache-level>=1 or 2. Alternatively one can set SM.caching=1,2 after the initialisation. In the command line mode, caching=n is equivalent to the number of "c" options (i.e., anyBSM SM -cc corresponds to <cache-level>=2 in the example above). Setting caching≥1 writes the found particle insertions for all generic topologies for a specific model to disk. The insertions for all Feynman diagrams are determined using a bruteforce algorithm. Therefore, this step can take up to several minutes for models with many particles. However, if caching≥1, this step is only required once. The option caching=2 (default) additionally caches the obtained analytic results to disk, such that the insertion of Feynman rules (i.e., definitions in the vertices.py and couplings.py) also have to be done only once. The analytic result saved to the cache is always expressed in terms of the abbreviations for the couplings defined in the UFO model files rather than the Feynman rules themselves.
The exact behaviour of the program depends on the chosen evaluation mode. If caching=2 but evaluation_mode='numerical' ('analytical'), the code internally first calculates the results using evaluation_mode='abbreviation' (if no result was found in the cache, otherwise the cached result is read from disk), saves it to the disk-cache and then evaluates it numerically (analytically). However, if caching<2 the individual diagrams are evaluated directly using the representations for the couplings/masses that correspond to the active evaluation mode.
One useful example choice may be if the program encounters an infra-red (IR) divergent loop function and returns a corresponding warning message. In some cases, this may not be critical if the loop function is multiplied by a vanishing coupling (i.e., the diagram does actually not exist but was computed because the information that the couplings vanish is not manifest in the UFO model). Thus, a good cross-check is to run the parameter point again without caching, which should yield the same numerical result but without a warning (since anyBSM does not calculate a diagram if any of the couplings vanish numerically).
F pyCollier pyCollier is a Python interface for the COLLIER Fortran library [77]. In the current pyCollier version, many but not all COLLIER functions are available. The pyCollier source code is hosted at https://gitlab.com/anybsm/pycollier.
Running the code requires at least Python version 3.5. The code is most easily used by installing the corresponding Python package by running pip install pyCollier which will automatically download and install pyCollier as well as most necessary dependencies. One necessary requirement which is not automatically handled by pip is the presence of a Fortran compiler (required for the compilation of the COLLIER library) such as gfortran or CLANG which can be installed from the systems system's package repository.
The pyCollier module is loaded via import pyCollier Loop integrals can then be evaluated e.g. via pyCollier.set_renscale(125**2) pyCollier.a0(125**2) where in the first step the renormalisation scale is set to 125 GeV. In the second step, the value of the A 0 scalar integral for an internal mass of 125 GeV is calculated. A detailed documentation of all available functions can be found at https://anybsm.gitlab.io/pycollier/pyCollier.html.