Two-loop Prediction of the Anomalous Magnetic Moment of the Muon in the Two-Higgs Doublet Model with GM2Calc 2

We present an extension of the GM2Calc software to calculate the muon anomalous magnetic moment ($a_\mu^{\text{BSM}}$) in the Two-Higgs Doublet Model. The Two-Higgs Doublet Model is one of the simplest and most popular extensions of the Standard Model. It is one of the few single field extensions that can give large contributions to $a_\mu^{\text{BSM}}$. It is essential to include two-loop corrections to explain the long standing discrepancy between the Standard Model prediction and the experimental measurement in the Two-Higgs Doublet Model. The new version GM2Calc 2 implements the state of the art two-loop calculation for the general, flavour violating Two-Higgs Doublet Model as well as for the flavour aligned Two-Higgs Doublet Model and the type I, II, X and Y flavour conserving variants. Input parameters can be provided in either the gauge basis or the mass basis, and we provide an easy to use SLHA-like command-line interface to specify these. Using this interface users may also select between Two-Higgs Doublet Model types and choose which contributions to apply. In addition, GM2Calc 2 also provides interfaces in C++, C, Python and Mathematica, to make it easy to interface with other codes.


Introduction
The anomalous magnetic moment of the muon, a µ , is one of the most important observables in particle physics. Its importance stems from the fact that it is one of the most precisely measured physical quantities and that it is very sensitive to new physics as well as the strong, weak and electromagnetic interactions of the Standard Model (SM). The significance of a µ has increased further in the wake of the recent release of the first results from the Fermilab Muon g − 2 experiment [1]. This collaboration found, in combination with the results from the Brookhaven National Laboratory [2], an experimental value of a Exp µ = (11659206.1 ± 4.1) × 10 −10 .
Phenomenological investigations of BSM physics are greatly enhanced by the use of precise software tools. These tools can automate the calculation of precision corrections, or provide numerical methods to quickly and reliably solve related problems. For example in the 2HDM one is often concerned that a particular benchmark point respects measured limits of electroweak oblique parameters, or that the new couplings and bosons do not provide contributions to Higgs decays which violate collider constraints. For the 2HDM a widely used software package is 2HDMC [113], which provides calculations of the spectrum, decays, oblique S, T and U parameters as well as a calculation of the 2HDM new physics contributions to a µ . There also exists the tool ScannerS [114] which can place theoretical, experimental, dark matter, and electroweak phase transition constraints on extended scalar sectors, including the 2HDM. A more precise calculation of the decays is available from the 2HDM decay dedicated code 2HDECAY [115], while PROPHECY4F [116] provides a package which focuses on calculating h → W W/ZZ → 4f decays. Additionally, the tools HiggsBounds [117][118][119][120][121] and HIGGSSIGNALS [122][123][124] allow one to place constraints on the Higgs sectors of BSM physics through the measured behaviour of Higgs bosons from collider search experiments.
Higher precision calculations for the spectrum, decays and other observables in the 2HDM are also available for arbitrary user-defined extensions of the SM through codes such as SARAH/SPheno [125][126][127][128][129][130] and FlexibleSUSY [131][132][133][134] (with decays recently added in FlexibleDecay [134]) where model files for the 2HDM are already distributed. Additionally, both packages provide one-loop contributions to a BSM µ . For two-loop contributions to a BSM µ in the MSSM, FlexibleSUSY links to the dedicated tool GM2Calc [135]. For the 2HDM there was no such option until now. Here we extend GM2Calc with the 2HDM to provide a program which includes state of the art two-loop level contributions of Ref. [102] to the anomalous magnetic moment of the muon.
The 2HDM contributions implemented in GM2Calc version 2 include the one-loop contributions, the two-loop fermionic corrections (including the well-known Barr-Zee diagrams), and the complete set of bosonic two-loop corrections. The one-loop contributions are of the order O(m 4 µ ) and therefore usually subdominant. The two-loop contributions arise at O(m 2 µ ) and are implemented at this order. Ref. [102] has assumed the flavour-aligned 2HDM with vanishing couplings λ 6,7 , but here we relax this assumption and allow the more general case. The 2HDM version of GM2Calc 2 allows the user to input deviations from the aligned limit, as well as select any one of the 2HDM types I, II, X or Y. Just like version 1, GM2Calc 2 allows the user to input parameter information using an SLHA-like [136,137] input file. We also provide interfaces in C, C++, Python and Mathematica to make it easy to link to other public codes. Furthermore, GM2Calc 2 (hereafter GM2Calc) can be used as a standalone tool for studies of a µ in the 2HDM, or to explore the 2HDM phenomenology more broadly it can be used in combination with other codes via the SLHA interface. For example, GM2Calc can be called alongside FlexibleSUSY, or in combination with other standalone tools like 2HDECAY or alongside the 2HDMC package, replacing its native calculation of the anomalous magnetic moment of the muon. The MSSM calculation is already available in GAMBIT [138,139] and it should be straightforward to extend this to also use the 2HDM calculation in future versions.
The rest of this paper is divided into the following sections. Section 2 provides a pedagogical introduction to the 2HDM, giving the Lagrangian both before and after electroweak symmetry breaking, and defining the Yukawa couplings in the various types of the 2HDM. There we give the one-and twoloop contributions to a BSM µ , and the uncertainty estimate for the calculation. Section 3 describes various ways to use the 2HDM in GM2Calc. It also shows how to interface GM2Calc using C++, C, Python and Mathematica. Finally, Section 4 demonstrates various possible ways GM2Calc can be applied to calculate the BSM contributions to the anomalous magnetic moment of the muon in the 2HDM.

The Model
The Two-Higgs Doublet Model (2HDM) extends the Standard Model (SM) with an additional scalar SU(2) doublet. We denote the two complex Higgs SU(2) doublets in the 2HDM as Φ i (i = 1, 2), where b i and c i are real scalar fields and a + i are complex scalar fields. Each Higgs doublet acquires a real non-zero vacuum expectation value (VEV), v i , which satisfy with 0 ≤ β ≤ π/2 and v ≈ 246 GeV, which implies v 1 = v cos β and v 2 = v sin β. The extended Lagrangian includes the Higgs potential L Scalar and the Yukawa interaction part L Yuk , The most general form of Higgs potential is where λ 5 , λ 6 , λ 7 and m 2 12 are real for the CP -conserving Higgs potential. The mass eigenstates of Higgs and Goldstone bosons h, H, A, H + , G 0 and G + are obtained as with from Eq. (7). The orthogonal matrices Z H 0 , Z A and Z H ± for Eq. (8) are given as The corresponding mass square matrices M 2 H 0 , M 2 A and M 2 H + for H 0 , A and H + in Eq. (9), respectively, are diagonalized as The Higgs boson mass eigenvalues are denoted as where in Feynman gauge m G + = m W and m G 0 = m Z . We chose the CP -even Higgs boson mixing angle α such that −π/2 ≤ β − α ≤ π/2 [113]. The general form of the Yukawa interaction Lagrangian is given as , and e 0 R are fermion gauge eigenstates. 2 In general the Yukawa coupling matrices, Γ 0 f and Π 0 f (f = u, d, l) are complex and non-diagonal 3 × 3 matrices which can not be simultaneously diagonalized, which produces flavour changing neutral currents (FCNC).
From the Yukawa interaction Lagrangian in Eq. (13) the 3×3 fermion mass matrices M f (f = u, d, l) are obtained as The fermion mass matrices are diagonalized using singular value decomposition. The diagonal fermion mass matrices M D f (f = u, d, l) are given by where V f and U f are unitary matrices, and the corresponding fermion mass eigenstates f = u, d, l are obtained as By combining Eq. (5) and Eqs. (14)- (16), Γ 0 u , Γ 0 d and Γ 0 l can be expressed in terms of M f and Π 0 f as The Yukawa interaction Lagrangian Eq. (13) can be also written in the so-called Higgs basis, where only one Higgs doublet has non-zero VEV. This can be obtained by rotating Φ 1,2 and the Yukawa matrices as such that only Φ v has the VEV v and Σ 0 f are its Yukawa couplings. It contains the SM-like Goldstone bosons, whereas the second field Φ ⊥ contains non-SM Higgs bosons H ± and A, and its Yukawa couplings ρ 0 f are important parameters in the 2HDM. 3 With this transformation the Yukawa interaction Lagrangian becomes (24) In this way, the fermion masses are purely generated from the Higgs doublet Φ v , and the fermion mass matrices in Eqs. (14)- (16) become • the 2HDM with general Yukawa structures.
We will now describe these variants of the 2HDM. The most common way to avoid FCNC is to impose a Z 2 symmetry, which leads to four different types of Yukawa interactions in which specific subsets of Yukawa couplings vanish. In the type I model all quarks and charged leptons couple to Φ 2 , so we set In the type II model all up-type quarks couple to Φ 2 , while the down-type quarks and charged leptons couple to Φ 1 , so we set In the type X model all quarks couple to Φ 2 , while the charged leptons couple to Φ 1 , so we set Type I Type II Type X Type Y In the type Y model the up-type quarks and charged leptons couple to Φ 2 , while the down-type quarks couple to Φ 1 , so we set In all these cases, trivially Γ 0 f and Π 0 f can be diagonalized simultaneously, and Higgs-mediated FCNC is avoided.
In the flavour-aligned Two-Higgs Doublet Model (FA2HDM) [83,84] it is assumed that the Π 0 f matrices are proportional to the Γ 0 f matrices, and we set where ζ f are the alignment parameters which are constrained by experimental results and used in phenomenological studies. Note that the aligned model contains the type I, II, X and Y models as special cases when the alignment parameters ζ f take the values given in Table 1.
We can summarize the fundamental Yukawa coupling matrices ρ f (f = u, d, l) for the different 2HDM variants and in different parametrizations as follows: for type I, II, Y, Y (using Table 1) and for the exact FA2HDM, Here the ζ f are the parameters introduced in Eq. (32) and Table 1, and the ( * ) notation indicates the conjugate is present for the case f = u and not present for f = d, l. GM2Calc offers two parametrizations for the general 2HDM. The "FA2HDM parametrization" starts from the FA2HDM and allows to directly modify the ρ f by additional matrices ∆ f , which represents the deviation from the flavour-aligned (or type I, II, X, Y) limit. The second "Π parametrization" starts from Eq. (27) and allows to directly specify the fundamental Yukawa matrices Π f . After suitable unitary transformations, the Yukawa interaction in the Lagrangian of Eq. (13) takes the following form: with the CKM matrix V CKM = V u V † d and the fermion mass eigenstates f = f L + f R , (f = u, d, l), defined in Eq. (20). The coupling matrices y S f are defined as [102]: and the coupling between the charged Higgs and each of the fermions is: where ρ f is given in Eq. (33).

Implemented contributions to a µ
Relevant 2HDM contributions to a µ arise at the one-loop and the two-loop level. The one-loop contributions are suppressed by two additional powers of the muon Yukawa coupling and thus of O(m 4 µ ); they are typically subleading unless some of the non-SM-like Higgs bosons have very small mass around a few GeV. The two-loop contributions can be classified into fermionic and bosonic contributions. They arise at the O(m 2 µ ) and are typically dominant. An important subclass of twoloop contributions are the so-called Barr-Zee diagrams, which contain a γγ-Higgs one-loop subdiagram. These diagrams have been studied extensively [140][141][142][143][144] and fully calculated in Ref. [100]. The full set of two-loop bosonic and fermionic diagrams (at O(m 2 µ )) has been computed in Ref. [102] (see also Ref. [103] for further phenomenological discussions of the individual contributions).
GM2Calc implements the full set of one-loop BSM contributions a 1 µ of O(m 4 µ ) and the two-loop BSM contributions a 2 µ of O(m 2 µ ) in the 2HDM from Ref. [102]. The full BSM contribution in the 2HDM (i.e. the difference between the 2HDM and the SM contributions) calculated by GM2Calc is therefore the sum In the following subsections we provide the explicit analytic results implemented in GM2Calc.

One-loop contributions
The one-loop BSM contributions to a 1 µ in the 2HDM is of O(m 4 µ ) and are given by with m l = (m e , m µ , m τ ), m ν = (m νe , m νµ , m ντ ) and Note, that the one-loop contribution from the SM Higgs boson h SM is subtracted in Eq. (40) and is therefore not included in a 1 µ . The loop functions F C 1 , F C 2 and F N 1 are given by

Two-loop contributions
The two-loop contributions of O(m 2 µ ) are divided into fermionic and bosonic loop contributions according to Ref. [102] as where the very small contribution from the shift of the Fermi constant, a ∆r-shift µ , is neglected. We generalize the fermionic two-loop contributions from Ref. [102] in the general 2HDM to the case of CKM mixing as follows: Note, that the two-loop contribution from the SM Higgs boson h SM is subtracted in Eq. (51) and is therefore not included in a 2 µ . The loop functions f S f (S = h, H, A, h SM ) are defined as: where is defined as: where f j is the SU(2) L partner of generation j of the fermion f i from generation i, and and The The bosonic two-loop contributions to the general 2HDM are composed of three parts as follows: The two terms a B,EW add µ and a B,nonYuk µ correspond to diagrams with the SM-like Higgs boson and diagrams without the new Yukawa couplings, respectively. They are typically subdominant, and full details on their definition and phenomenological impact are given in Refs. [102] and [103]. The Yukawa contribution is given by where the equations for a η λ,t (with the common prefactor put in the front of the above Eq. (67)) are given in the appendix of Ref. [102]. Compared to this reference, the expression has been generalized to include λ 6 and λ 7 . 4 These bosonic contributions are implemented in the realistic approximation of small cos(β − α), and only terms up to linear order in cos(β − α) are taken into account. In this scenario the boson h has the tree-level couplings of the SM Higgs boson [146]. Furthermore, the bosonic contributions are derived only for the flavour-aligned 2HDM. Referring to the different cases of Yukawa couplings in Eq. (33), the bosonic corrections apply for the cases of the discrete symmetries, for the FA2HDM, and also for the general 2HDM in FA2HDM parametrization (assuming the ∆ f matrices are small). In contrast, we set ζ l = 0 in Eq. (67) in the case of the general 2HDM in Π parametrization since in this case the bosonic corrections are not applicable in this form.
The bosonic two-loop contributions are also available in the GM2Calc source code and in the form of a Mathematica file.

Running couplings
Among the fermionic two-loop contributions, the diagrams with an internal top quark, bottom quark or tau lepton loop give the dominant contribution to a BSM µ . These diagrams are proportional to the values of the fermion masses in the loop. So long as there is no three-loop calculation available in the 2HDM, it is formally irrelevant which renormalization scheme to use for the fermion masses in the This structure clarifies the appearance of Λ 567 in Eq. (67).
loop. In GM2Calc two possible definitions of these fermion masses are implemented: Running masses: In the "input masses" scheme, the top quark pole mass m t , the MS bottom quark mass in the SM with five active quark flavours, m MS b , at the renormalization scale Q = m MS b , and the tau lepton pole mass m τ are used in the loops. These masses are typically used as input for spectrum generators [136,137]. In the "running masses" scheme, the running MS top quark mass m MS t (Q), the MS bottom quark mass m MS b (Q) and the MS tau lepton mass m MS τ (Q) are used in the loops. The renormalization scale Q is set to the mass of the Higgs boson in the Feynman diagram. The "running masses" scheme is also used in 2HDMC [113]. The difference between these two schemes is shown in Section 4.

Uncertainty estimate
We provide an uncertainty estimate for a 2 µ as follows: where and The term δa 2 ,∆r µ accounts for the fact that the two-loop contribution a ∆r-shift µ has been neglected in Eq. (49). In Refs. [102,147] it was shown that in the relevant parameter space |a ∆r-shift µ | ≤ 2 × 10 −12 , which we use as an upper bound in the uncertainty estimate. The term δa 2 ,m 4 µ µ estimates missing two-loop terms of O(m 4 µ ) using the known universal two-loop QED logarithmic contributions [148]. The term a 3 µ is an estimate for the expected three-loop contributions. From experience the QED contributions are among the largest contributions at each loop level. For this reason we use again the known universal logarithmic QED contributions from Ref. [148] to provide an estimate for the unknown three-loop contributions.

Requirements
To build GM2Calc the following programs and libraries are required:

Running GM2Calc from the command line
GM2Calc can be run from the command line with an SLHA-like [136,137] input file as where example.thdm is the name of the SLHA-like input file. Alternatively, the input parameters can be piped in an SLHA-like format into GM2Calc as The calculated value for a BSM µ and δa µ is written to the standard output. Depending on the input options and parameters, the output may contain the following block with a BSM µ and δa µ : The SLHA-like format for the input parameters is defined in the following subsections. 6 For instructions and examples of running the MSSM evaluation see Ref. [135].

General options
In the SLHA-like input the selection of the configuration flags for GM2Calc can be given in the GM2CalcConfig block as defined in Ref. [135]. An example GM2CalcConfig block reads as follows: # calculate uncertainty (0 or 1)

Standard Model input parameters
The Standard Model input parameters are read from the SMINPUTS, GM2CalcInput and VCKMIN blocks. Example blocks that define the Standard Model input parameters may read:   The entries of the SMINPUTS and VCKMIN are defined in Refs. [136,137]. In the block entry GM2CalcInput [33] the mass of the Standard Model Higgs boson can be given. Unset parameters in the SMINPUTS and GM2CalcInput blocks are assigned default values. Unset parameters in the VCKMIN block are assumed to be zero.

Two-Higgs Doublet Model input parameters
We allow the specification of the 2HDM input parameters in two different "bases", see Table 2. The basis parameters are defined as in [113]. In the "gauge basis" the Lagrangian parameters, λ 1,...,7 are used as input. In the "mass basis" the Higgs boson masses and the mixing parameter sin(β − α) are used as input, instead of the Lagrangian parameters λ 1,...,5 , but λ 6 and λ 7 can still be used. All available input parameters are listed in Table 3.
Mass basis input parameters. The "mass basis" input parameters are read from the MINPAR and MASS blocks for compatibility with 2HDMC [113]. The following shows an example input in the "mass basis" for the type II 2HDM: MHp R ≥0 parameters common to both the gauge and mass basis Unset parameters in the MINPAR and MASS blocks are assumed to be zero, except for tan β which will raise an error. Entry 24 of the MINPAR is used to select the type of 2HDM. Specifically the integer values 1, . . . , 6 for the Yukawa type correspond to: 1 = type I, 2 = type II, 3 = type X, 4 = type Y, 5 = Flavour-Aligned, 6 = general 2HDM. The input entries 21, 22, and 23 in the MINPAR block can be used to set the values of ζ u , ζ d and ζ l in the Flavour-Aligned 2HDM, and are ignored for all other types. Additional blocks for the general case are described below.
Gauge basis input parameters. Alternatively, the input can be given in the "gauge basis" in a format compatible with 2HDMC [113]. In the "gauge basis" the 2HDM parameters must be given in the MINPAR block. The available "gauge basis" input parameters are listed in Table 3. The following shows an example input in the "gauge basis": Two parametrizations for the general 2HDM are implemented, see Eq. (33). The first option is a deviation from the FA2HDM (or types I, II, X, Y), parametrized by the additional matrices ∆ f . Thus, after choosing Yukawa type = 1, . . . , 5, i.e. type I, II, X, Y and FA2HDM, the matrices ∆ u , ∆ d and ∆ l can be optionally given in the following dedicated blocks: The other option for the general 2HDM corresponds to the Π parametrization, implemented when setting Yukawa type = 6. In this case, the input parameters ∆ u , ∆ d and ∆ l are ignored and the real parts of the matrices Π u , Π d and Π l can be given in the following dedicated blocks: 1 Block G M 2 C a l c T H D M P i u I n p u t 2 The input parameters Π u , Π d and Π l are ignored when choosing Yukawa type = 1, . . . , 5, i.e. type I, II, X, Y and FA2HDM.
The SLHA-like interface described so far is very convenient and intuitive to use and it does not require any prior knowledge of programming languages to run GM2Calc this way. The execution time for this is also very short, around 5 ms per point on current laptops. 7 Nevertheless even faster execution times and easy interfacing with other existing calculators and sampling algorithms are enabled through our C++ (0.05 ms), C (0.05 ms), Mathematica (0.5 ms) and Python (0.08 ms) interfaces. In the following subsections we describe how to use each of these interfaces.

Running GM2Calc from within C++
GM2Calc provides a C++ programming interface, which allows for calculating of a BSM µ in the 2HDM up to the two-loop level. The following C++ source code snippet shows a two-loop example calculation in the 2HDM of type II with input parameters defined in the "mass basis", see Table 2.
1 # include " gm2calc / gm2_1loop . hpp " 2 # include " gm2calc / gm2_2loop . hpp " 3 # include " gm2calc / g m 2_ un ce r ta in ty . hpp " 4 # include " gm2calc / gm2_error . hpp " 5 # include " gm2calc / THDM . hpp " 6 7 # include < cstdio > This example source code can be compiled as follows (assuming GM2Calc has been compiled to a shared library libgm2calc.so on a UNIX-like operating system with g++ installed): Here example.cpp is the file that contains the above listed source code. The variable GM2CALC_DIR contains the path to the GM2Calc root directory and EIGEN_DIR contains the path to the Eigen library header files. Running the created executable a.out yields 1 $ ./ a . out 2 amu = 1.67323 e -11 + -3.3616 e -12 In line 12 of the example source code an object of type Mass_basis is created, which contains the input parameters in the mass basis, see Table 2. The mass basis input parameters are set in lines 13-31. Note that the Yukawa type is set to type II in line 13, which implies that given values of ζ f are ignored and internally fixed to the values given in Table 1. Additionally the inputs Π f and ∆ f are unused for this type, and ignored. In line 34 an object of type SM is created that contains all SM input parameters. The SM input parameters are set to reasonable default values from the PDG [154]. Note that the 2HDM model should be created within a try block, because the constructor of the 2HDM class throws an exception if a physical problem occurs (e.g. a tachyon) or an input parameter has been set to an invalid value, see Table 3.
Alternatively, the 2HDM input parameters can be given in the "gauge basis", i.e. in terms of the Lagrangian parameters, see Table 2. The following C++ source code snippet shows a corresponding example two-loop calculation in the 2HDM of type II, where the input parameters are defined in the "gauge basis".

Running GM2Calc from within C
Alternatively to the C++ programming interface detailed in Section 3.4, GM2Calc also provides a C programming interface. The following C source code snippet shows a two-loop example calculation in the 2HDM of type II with input parameters defined in the "mass basis", see Table 2.
The C example source code is very similar to the mass basis C++ example. In line 12 an object of type gm2calc_THDM_mass_basis is created and filled with the mass basis input parameters in lines 13-32. The Yukawa type of the model is defined in line 13 to be type II. In line 35 an object of type gm2calc_SM is created, which contains the SM input parameters. The Alternatively, the 2HDM input parameters can be given in the "gauge basis", similarly to the gauge basis C++ example. The following C source code snippet shows a corresponding example two-loop calculation in the 2HDM of type II, where the input parameters are defined in the "gauge basis".

48 Print [ result ];
In line 1 GM2Calc's MathLink executable bin/gm2calc.mx, which is created when building GM2Calc, is loaded into the Mathematica session. In lines 4-6 two configuration options to customize the calculation are set: The calculation shall be performed at the two-loop level using the "running masses" scheme defined in Section 2.2.3. In lines 9-29 the SM input parameters are defined. Unset parameters are set to reasonable default values, see Options[GM2CalcSetSMParameters]. In lines 32-46 the values of a BSM µ and δa µ are calculated using the function GM2CalcAmuTHDMGaugeBasis, which takes the gauge basis input parameters as arguments, see Table 3. The result is printed in line 48.

Running GM2Calc from within Python
Newly implemented in GM2Calc 2.0.0 is the ability to interface with Python using the package cppyy. An example calculation using the interface is shown in the code snippet below, working in the "mass basis". Note similarity between the above code and the C++ and C interfaces. Line 3 is to make sure that this example which is written using Python 3-style print functions can still work in Python 2. Line 4 imports the interface script gm2_python_interface which loads the cppyy and os packages, as well as the header and library locations. This interface file is originally in the src subdirectory. After performing the interface script will be copied into the subdirectory bin, and it will be filled with the path information for the GM2Calc headers, library, and the Eigen3 path. The interface can be imported from there by other Python scripts, or moved to an appropriate location where the user has their own Python scripts. Lines 6-11 load the relevant header files, and line 13 loads the GM2Calc shared library. In lines 16-20 the necessary namespaces from C++ are loaded into Python. In line 23 the 2HDM Mass_basis object is initialized, while on line 24 the 2HDM is specified to be type II. Lines 25-36 involve setting the values for simple attributes in the basis. Lines 36-42 assign values to the Eigen::Matrix attributes, however since these are meant to be ignored, they are just set to 0. Lines 44-49 initialize an SM object and ensures it has the appropriate parameters. Line 51 initializes the config object, which is used to flag the use of running coupling in the next line. Line 56 initializes a THDM object using the Mass_basis, SM, and Config information. Lines 58-63 prints out the values of a BSM µ and δa µ which are calculated using the interface functions calculate_amu_1loop, calculate_amu_2loop, and calculate_uncertainty_amu_2loop. Alternatively an error message will be printed out on line 63 should a problem arise.
Another example of the Python interface is shown below, this time using the "gauge basis": 1 # !/ usr / bin / env python 2 3 from __future__ import print_ function 4 from g m 2 _ p y t h o n _ i n t e r f a c e import * 5 6 cppyy . include ( os . path . join ( " gm2calc " ," gm2_1loop . hpp " ) ) 7 cppyy . include ( os . path . join ( " gm2calc " ," gm2_2loop . hpp " ) ) In line 23 we instead initialize a Gauge_basis object. To define the attribute lambda, we need to circumvent Python's reserved keywords. This is done by defining a 7 ×

Parameter scan in the type II and X models
As an application we perform a 2-dimensional parameter scan over m A and tan β for the type II and type X 2HDM models, similarly to Ref. [85]. However, in contrast to Ref. [85] we include the two-loop bosonic contributions and use the updated value of a BSM µ = (25.1 ± 5.9) × 10 −10 from Eq.
(3). The following C++ source code shows the program to perform the scan.

Size of fermionic and bosonic contributions
In the following we illustrate the calculation of the two-loop fermionic and bosonic contributions, a F µ and a B µ , separately. For the illustration we perform a scan over m A for the demo parameter scenario from 2HDMC [113], which is a type II 2HDM scenario where m H = 400 GeV, m H ± = 440 GeV, tan β = 3, sin(β − α) = 0.999, λ 6 = λ 7 = 0 and m 2 12 = (200 GeV) 2 . The following C source code shows the program to perform the scan. . yukawa_type = gm2calc_THDM_type_2 , 18 . mh = 125 , 19 . mH = 400 , 20 . mA = mA_start + i *( mA_stop -mA_start ) / N_steps , 21 . mHp = 440 , 22 . s i n _ b e t a _ m i n u s _ a l p h a = 0.999 , 23 . lambda_6 = 0 , 24 . lambda_7 = 0 , 25 . tan_beta = 3 , 26 . m122 = 40000 , 27 . zeta_u = 0 , 28 . zeta_d = 0 , 29 . zeta_l = 0 ,  [103]. It should be noticed that, once the masses for the charged and CP-even scalar are close together and below around 300 GeV, electroweak as well as unitarity and perturbativity constraints can be evaded for arbitrarily low values of m A [85]. The input parameters here are chosen accordingly. Since m A < m h SM /2, it is possible for h → AA decays to occur unless we enforce the coupling C hAA = 0. 8 This fixes the value of λ 1 according to Eq. (12) in Ref. [103], which can be set in the mass basis using m 2 12 and applying the relations in Eqs. (2.12)-(2.13) in Ref. [102]. This leads to the fitted 2nd-order polynomial relation with a dependence on m A seen in the source code below.
These scenarios have a very light pseudoscalar mass, but LHC limits are much weaker compared to the type II case and can be evaded for these scenarios. The two-loop fermion contributions rise rapidly as the pseudoscalar mass decreases, dominating over the two-loop bosonic contributions, though the latter are just large enough to have an impact on constraints from a µ . Note that for higher values of m A it is possible to get larger bosonic contributions as can be seen in Figure 10 of Ref. [103]. In the scenarios we plot here the one-loop contributions, which are not shown in Figure 2, have a negative effect on the contributions, with a size of approximately one-third of the two-loop fermionic contributions. Thus it can be seen that a BSM µ can be explained with a very low m A for the values of m A in the purple region. The scan for this scenario can be performed with the following C source code: 1 # include " gm2calc / gm2_1loop . h " 2 # include " gm2calc / gm2_2loop . h " 3 # include " gm2calc / gm2_error . h " 4 # include " gm2calc / THDM . h " 5 6 # include < stdio .h > 7 # include < math .h > 8 9 int main () 8 The coupling C hAA has been discussed in detail in Ref. [103]. Setting C hAA = 0 directly corresponds to the choice Λ 6 = 0, where Λ 6 is a potential parameter in the so-called Higgs basis, and to a specific relation between all potential parameters in the general basis. It has been checked that although non-null values for C hAA can be allowed, they are experimentally strongly constrained. Thus, for simplicity, we have adopted C hAA = 0 in our analysis, which automatically guarantees that λ 1 (or equivalently m 2 12 ) will not be chosen in a region excluded by experimental constraints or by unitarity or perturbativity.   . yukawa_type = gm2calc_THDM_aligned , 21 . mh = 125 , 22 . mH = 150 , 23 . mA = mA , 24 . mHp = 150 , 25 . s i n _ b e t a _ m i n u s _ a l p h a = 0.999 , 26 . lambda_6 = 0 , 27 . lambda_7 = 0 , 28 . tan_beta = 2 , 29 // Polynomial fit following Eq . (2.12) from PRD 30 . m122 = 3187.3 + mA *(3.27803 + 0.0165557* mA ) , 31 . zeta_u = -0.1 , 32 . zeta_d = -0.1 , 33 . zeta_l = 50 , 34 .

Running fermion masses
In this subsection we study the effect of using the input vs. running fermion masses in the twoloop fermionic contributions as described in Section 2.2.3, and compare the results with 2HDMC. The following Mathematica source code shows a program to perform a scan over m A using the same type II 2HDM parameter region as in Section 4.2. The function CalcAmu calculates the two-loop fermionic contribution a F µ and the uncertainty δa µ for a given value of m A and the mass basis input parameters defined above. In line 20 the usage of running fermion masses is disabled and the calculation is performed in the subsequent line. Similarly, in line 24 the usage of running fermion masses is enabled and the calculation is performed in the subsequent line. The results are collected in the variable data and are exported to a file in line 31. We also used a very similar script to do the same calculations for the FA2HDM scenarios discussed in Section 4.2.
The effect of using running fermion masses in the two-loop fermionic contributions is shown in Figure 2 for the type II (bottom left panel) and flavour aligned (bottom right panel) scenarios matching those described in the previous section. The red solid line shows the value of a F µ when the running masses (73) are used, i.e. the 3rd generation fermion masses are run to the scale of the Higgs boson in the two-loop fermionic Barr-Zee Feynman diagrams. Note that although the vertical axes are slightly different, the red lines shown in the bottom panels are identical to the red lines from the corresponding panels immediately above them, which were discussed in the previous section. The green dashed lines show the value of a F µ when input fermion masses listed in (72) are used. In both scenarios that we look at the value of a F µ is smaller when running masses are used, though the difference is only distinguishable for the type II case where the size of the contributions is much smaller. The reason for this is that due to the negative fermion mass β functions the running masses are numerically smaller than the corresponding input masses in the shown scenarios, which leads to a systematic reduction of the fermionic two-loop contributions.
In addition to the red solid and green dashed lines from GM2Calc, as the scenario in question is a benchmark point for 2HDMC [113], we show in the bottom left panel of Figure 2 the corresponding result obtained with 2HDMC 1.8.0 as black dotted line. Note that at the two-loop level 2HDMC includes only fermionic contributions to a BSM µ . Furthermore, 2HDMC does not subtract the contributions from the SM Higgs boson and does not include the two-loop contributions from the charged Higgs boson. Therefore to obtain this black dotted line we have thus subtracted the two-loop SM Higgs contributions from the 2HDMC result and added the two-loop contributions from the charged Higgs boson. Since 2HDMC inserts running fermion masses into the fermionic contributions, the black dotted line can be compared to the red solid line in the figure. There is a small deviation between these two lines, which originates from the inclusion of fermionic Barr-Zee diagrams with an internal Z boson in GM2Calc, which are not included in 2HDMC.
In the bottom panels of Figure 2 we also show the uncertainties calculated with (74) as lighter shaded regions of the corresponding color about the red and green lines. In the bottom left panel the red and green lines both lie within the uncertainty estimate for the alternative prediction (shaded green and shaded red regions respectively) for all values of m A plotted. This is also true in the bottom right panel though there is no visible distinction between the red and green lines or their uncertainties here. Since the difference between the lines is of higher order, this indicates that our uncertainty estimate is working as expected and accounts for the expected higher order corrections.

Summary
We have presented version 2 of GM2Calc, with its new capability to calculate the BSM contributions to the anomalous magnetic moment of the muon in the 2HDM. The contributions include all the oneloop diagrams, two-loop fermionic Bar-Zee diagrams, as well as the bosonic two-loop contributions. The new version of GM2Calc provides the calculation of the 2HDM contributions with a precision of up to O(m 4 µ ) at one-loop and O(m 2 µ ) at two-loop level along with an estimate of uncertainty of the evaluation. GM2Calc performs this state of the art precision calculation at high speed, with execution times that can be as short as O(0.05 ms) per point, allowing for rapid sampling of the parameter space.
GM2Calc is easy to configure and run. The user can select well-known types of the 2HDM, specifically type I, II, X, Y as well as the flavour-aligned version (FA2HDM), or the fully general 2HDM. For the latter the user can specify the inputs as deviations away from the flavour alignment of the FA2HDM (or from type I, II, X, and Y) or by directly specifying the more fundamental Yukawa matrices, Π f defined in Eq. (27). The user can also decide whether they will give inputs in the gauge basis using λ 1,...,7 or the mass basis using m h,H,A,H ± , λ 6,7 , and sin(β − α). The input parameters and settings can be specified in an SLHA-like input file, mirroring the original MSSM version. Additionally, GM2Calc can be interfaced to other programs using C++, C, Mathematica, or Python, the latter being a new interface developed for GM2Calc 2.
For each of these interfaces we presented simple and easy to follow usage examples in Section 3, that are straightforward to adapt. In Section 4 we have also presented some applications and results, demonstrating different features of the code. In each of these we show the source code in the manual and also provide them as supplementary files. GM2Calc is actively developed on GitHub and users with any questions may contact the authors through our GitHub page or directly by email.