Likelihood analysis of the flavour anomalies and $g-2$ in the general two Higgs doublet model

We present a likelihood analysis of the general two Higgs doublet model, using the most important currently measured flavour observables, in view of the anomalies in charged current tree-level and neutral current one-loop rare decays of $B$ mesons in $b\to c l \overline{\nu}$ and $b\to s\mu^{+}\mu^{-}$ transitions, respectively. We corroborate that the model explains the latter and it is able to simultaneously fit the experimental values of the $R(D)$ ratio at $1\sigma$, but it can not accommodate the $R(D^{*})$ and $F_{L}(D^{*})$ observables. We find that the fitted values for the angular observables in $b\to s\mu^{+}\mu^{-}$ transitions exhibit better agreement with the general two Higgs double model in comparison to the SM. We also make predictions for future collider observables $\mathrm{BR}(t\to ch)$, $\mathrm{BR}(h\to bs)$, $\mathrm{BR}(h\to \tau\mu)$, $\mathrm{BR}(B_{s}\rightarrow\tau^{+}\tau^{-})$, $\mathrm{BR}(B^{+}\rightarrow K^{+}\tau^{+}\tau^{-})$ and the flavour violating decays of the $\tau$ lepton, $\mathrm{BR}(\tau\rightarrow3\mu)$ and $\mathrm{BR}(\tau\to\mu\gamma)$. The model predicts values of $\mathrm{BR}(t\to ch)$, $\mathrm{BR}(B_{s}\rightarrow\tau^{+}\tau^{-})$ and $\mathrm{BR}(B^{+}\rightarrow K^{+}\tau^{+}\tau^{-})$ that are out of reach of future experiments, but its predictions for $\mathrm{BR}(h\to bs)$ and $\mathrm{BR}(h\to \tau\mu)$ are within the future sensitivity of the HL-LHC or the ILC. We also find that the predictions for the $\tau\rightarrow3\mu$ and $\tau\to\mu\gamma$ decays are well within the projected limits of the Belle II experiment. Finally, using the latest measurement of the Fermilab Muon $g-2$ Collaboration, we performed a simultaneous fit to $\Delta a_{\mu}$ constrained by the charged anomalies, finding solutions at the $1\sigma$ level. Once the neutral anomalies are included, however, a simultaneous explanation is unfeasible.


Introduction
The Standard Model (SM) of particle physics contains three fermion families which acquire mass by means of the interaction with the Higgs boson. The two Higgs doublet model (2HDM) is one of the simplest ways to extend the Higgs sector, which is the least constrained sector of the Standard Model. Two Higgs doublets also appear in many more elaborate extensions of the SM that are based on fundamental principles, such as supersymmetry (see e.g. [1]), the Peccei-Quinn symmetry [2,3] or grand unified theories (see [4] for a recent review). Two Higgs doublet models are also motivated from electroweak baryogenesis studies, where it has been shown that contributions coming from the new physical Higgs bosons to the effective Higgs potential can strengthen the phase transition and in addition introduce new sources of charge-parity (CP) violation, from both fermion and scalar sectors [5][6][7][8][9][10][11][12][13]. As a result, the 2HDM is one of the most popular SM extensions and has been frequently used as a benchmark for phenomenological studies (see e.g. [14] for a review of 2HDM studies). Furthermore, the presence of another doublet can contribute to resolving anomalies in lepton flavour universality observables [15,16] and muon g-2 , while scenarios where the extra doublet is "inert" can also explain dark matter [47][48][49][50][51][52][53][54][55][56][57][58][59].
The new interactions between the SM fermions and the physical states arising from the introduction of a second Higgs doublet imply a richer phenomenology than the SM. This is further enhanced by the new free parameters and couplings in the general two Higgs doublet model (GTHDM), also known as type-III 2HDM [60]. Physical effects such as CP violation, scalar mixing and flavour changing transitions are expected [61], allowing for signatures to be observed in particle colliders. One of the most interesting experimental consequences of the flavour changing currents present in the GTHDM is lepton flavour universality (LFU) violation. Experimental measurements of LFU violation come from flavour changing charged currents (FCCCs), such as those in B meson decays, and flavour changing neutral currents (FCNCs), for instance in kaon decays. The observed deviations from the SM in the measurements of FCCCs (around 3.1σ from the SM [62]) and FCNCs (close to a combined 6σ deviation, see for example [63][64][65]), hint at the existence of new physics (NP) contributions and thus serve as a clear motivation for the study of NP models capable of explaining the anomalies.
It has indeed been shown that the GTHDM is able to explain the charged anomalies at 2σ [15,16,66,67]. Similar analyses for the neutral anomalies have also been presented previously [15,[68][69][70], finding solutions at the 2σ level and up to the 1σ level including right-handed neutrinos [70]. Nevertheless, the majority of these studies have only explored solutions in restricted regions of the parameter space, with a lack of discussion of the role of (marginally) statistically preferred regions, and often considering the b → sll observables from model independent global fits [15,70]. Statistically rigorous explorations of the parameter space of the model contrasted directly to experimental constraints have rarely been performed, and even those were focused exclusively on interactions in the quark sector [71].
Furthermore, the longstanding discrepancy between the experimentally measured and SM predicted values of the anomalous magnetic moment of the muon a µ has recently been brought back to the spotlight with the new measurement by the Muon g-2 experiment at Fermilab [72]. The latest experimental value, taking into account the measurements at both Brookhaven National Laboratory and Fermilab, is a Exp µ = 116592061 ± 41 × 10 −11 . Compared to the theoretical prediction in the SM from the recent Muon g − 2 Theory Initiative White Paper, a SM µ = 116591810 ± 43 × 10 −11 [73], building on the extensive work examining the various SM contributions in [74][75][76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][92][93], the measured value differs from the SM prediction by ∆a µ = 2.51 ± 59 × 10 −9 , corresponding to a discrepancy of 4.2σ. Models with a second Higgs doublet have been studied extensively in the literature as sources to explain this deviation . However, no simultaneous global fit of the flavour anomalies and a µ in the GTHDM has been attempted giving a proper statistical insight into the whole parameter space. Therefore, in this paper we present a frequentist inspired likelihood analysis for the GTHDM, simultaneously including the FCCC observables, both b → sµ + µ − transitions and the muon anomalous magnetic moment, along with other flavour observables. We perform a global fit of all constraints using the inference package GAMBIT, the Global And Modular Beyond-the-Standard-Model Inference Tool [94,95]. GAMBIT is a powerful software framework capable of performing statistical inference studies using constraints from collider [96], dark matter [97], flavour [98] and neutrino [99] physics, as well as cosmology [100]. It has already been used for detailed statistical analyses of a variety of beyond the Standard Model (BSM) models, including supersymmetry [101][102][103], scalar singlet DM [104][105][106][107][108], axion and axion-like particles [109,110], and neutrinos [99,111], as well as an initial analysis of the 2HDM [112]. Our work enhances the FlavBit [98] and PrecisionBit [113] modules of GAMBIT to support the GTHDM. We also make use of various external codes: SuperIso 4.1 [114][115][116][117] for computing flavour observables, the 2HDMC 1.8 package [118] for precision electroweak constraints, the HEPLike package [119] which provides likelihoods for the neutral anomaly related observables, and the differential evolution sampler Diver 1.0.4 [120].
The paper is organised as follows. In section 2 we present the Higgs and Yukawa sectors along the theoretical bounds for their parameters. In section 3 we define the effective Hamiltonian and the Wilson coefficients (WCs) for b → sµ + µ − transitions. Then, in section 4 we list the observables to be used in our scans. Following this, our results from the global fit and predictions for future experiments in colliders are discussed in section 5. Finally, we summarise our conclusions in section 6.

GTHDM
The GTHDM has been actively investigated in both its scalar and Yukawa sectors. These can be written in three different ways, namely in the generic, Higgs and physical bases, all of them related via basis transformations [121]. Particularly, with respect to the Yukawa sector, in the past theorists imposed discrete symmetries to avoid flavour changing transitions, the most popular being the Z 2 symmetry in the type-II 2HDM [122,123]. However, it has been shown that there is no fundamental reason for forbidding flavour changing couplings [124]: if the mixing angle is small, the non-observation of several tree level flavour changing transitions can be explained by the alignment phenomenon. This, and a suppression inversely proportional to the mass of the heavy Higgses in the tree level amplitudes, could suppress the effects coming from the off-diagonal Yukawa couplings, without invoking the so called natural flavour conservation (NFC) condition [122].
Here we review the Higgs potential and the Yukawa Lagrangian of the model as well as the relevant theoretical constraints coming from stability, unitarity and perturbativity at leading order (LO). We also make use of the precision electroweak constraints from the oblique parameters. For a more comprehensive review of the model the reader is referred to [14,[125][126][127].

Higgs potential
The most general renormalizable scalar potential in the GTHDM is commonly written as [14,128] where the two scalar doublets are given by with υ i the vacuum expectation values (VEV) of the fields, while linear combinations of the fields ρ i , η i and φ ± i form mass eigenstates: where the fields φ + i are charged complex scalars. From the eight degrees of freedom, three of them (G W ± and G Z ) get absorbed by the longitudinal components of the vector bosons. The remaining five make up the new particle spectrum of the model, namely, h and H are physical CP-even states, A is a CP-odd state and H ± are two charged Higgs bosons. The rotation matrices are defined according to In this work, we assume a CP conserving scalar sector, which implies all the parameters in Eq. (2.1) to be real [128]. Additionally, for simplicity, we set λ 6 = λ 7 = 0. In particular, for this choice of the quartic couplings, the necessary and sufficient conditions to ensure positivity of the potential along all directions are given by [14,128] whereas the tree level unitarity of the couplings imposes [14] |a ± | , |b ± | , |c ± | , |d ± | , |e ± | , |f ± | < 8π, Following [71,125] we also include the oblique parameters S, T and U , which parametrise radiative corrections to electroweak gauge boson propagators. In this study we computed these oblique parameters with the 2HDMC package and these are contrasted with the most probable values inferred from experimental data, as found by the Gfitter group [129] S = 0.05 ± 0.11, T = 0.09 ± 0.13, U = 0.01 ± 0.11, (2.14) with correlations given by (2.15)

Yukawa Lagrangian
The most general Yukawa Lagrangian in the generic scalar basis {Φ 1 , Φ 2 } reads [71]: 16) where the superscript "0" notation refers to the flavour eigenstates, andΦ j = iσ 2 Φ † j . The fermion mass matrices are given by Notice that this matrices need to be diagonalized. This can be done through a bi-unitary transformationM where the fact that M f is Hermitian implies that V f L = V f R , and the mass eigenstates for the fermions are given by Then, Eq. (2.17) takes the formM Using the expressions above we can write the Yukawa Lagrangian in the mass basis as 1 where a, b = 1, 2, 3 and (2.26) At first, the total number of new complex Yukawa couplings to consider is 54. Considering only their real parts and the ansatz we get only 12 Yukawa parameters (i.e., ignoring 3 → 1 and 2 → 1 generation transitions).
Here, the ξ u matrix has been previously considered to be asymmetric from B s − B s oscillations constraints at one loop level and for heavy Higgs masses of order 700 GeV [130,131]. However, since we are approaching the dominant contribution process at LO and we are exploring masses in the range [0.5, 4.0] TeV as in [71], we will consider only the symmetric case, i.e., ξ u tc = ξ u ct . Hence, assuming the remaining ξ d and ξ l matrices to be symmetric as well, the total number of parameters to scan over is reduced by 3.

Effective Hamiltonians for flavour changing transitions
Most of the relevant flavour observables that we consider in this work arise from processes with either suppressed or negligible contributions from SM particles. Hence, these processes are often dominated by BSM contributions, which can be generated by a large variety of UV complete theories. It is often convenient to study these transitions using the, modelagnostic, effective Hamiltonian approach, where transition operators are decomposed using the Operator Product Expansion (OPE) into a collection of simple, low-energy, operators. Associated with each of these operators comes a WC, which encodes the knowledge of the high-energy theory. Even for complete high-energy theories, as it is our case, it is extremely useful to work with the effective Hamiltonian, since one can easily compute most observables of interest in terms of a small set of WCs. In fact, there are only two independent flavour changing transitions that give rise to the majority of the studied observables, and these are the neutral b → s + − transition and the charged b → c ν transition. In this section we write down the effective Hamiltonian for both of these transitions and provide expressions for the BSM contributions to the WCs that arise in our model.
where µ is the energy scale at which the WCs are defined, and are the FCNC local operators encoding the low-energy description of the high energy physics that has been integrated out. The prime operators are obtained by the replacement P R(L) → P L(R) . The WCs can be written as where C SM i is the SM contribution to the ith WC and ∆C i is the NP contribution, a prediction of the GTHDM model. The SM contribution to the scalar WCs, C as computed with SuperIso. We evaluate the NP scalar and pseudoscalar coefficients ∆C ( ) S,P at tree level, which is the LO contribution from the GTHDM [70]. Henceforth we will use the scalar and pseudoscalar coefficients in the basis defined in SuperIso, i.e., C The remaining coefficients, ∆C 7,8,9,10 first appear at one loop level and we therefore include the one-loop BSM contributions to these in our analysis. These one-loop corrections can be split by contribution as follows,

7)
∆C 9 = C γ 9 + C Z 9 + C box 9 , (3.8) where C Z 9,10 and C γ 7,9 come from the Z and γ penguins, respectively (figure 1), and C box 9,10 are contributions from box diagrams, (figure 2). At this level, the ∆C 9 and ∆C 10 coefficients are suppressed as m b /m t with respect to their non-prime counterparts. However, for studying the effects of flavour-changing Yukawa couplings we include these coefficients for completeness. C g 8 is the WC related to the chromomagnetic operator coming from gluon penguins and the NP contributions ∆C 7,8 are computed in [70].

Penguins and boxes computation
We review the computation of the WCs in Eqs. (3.7-3.9) which have been obtained already for both the flavour conserving general THDM in SuperIso and for the GTHDM itself in [15,70,132]. In these latter works, the Yukawa couplings related to ξ d were assumed to be zero or negligibly small from the beginning, avoiding the appearance of possible mixed terms between the down and up couplings that, at first, might not be as suppressed as those involving only down quarks. This computation is performed assuming = µ in the final state, as inspired by our choice of Yukawa textures in Eq. (2.27), but it can be easily generalised for all flavours when required.
Using the model files provided by FeynRules from [133], we generate in FeynArts the one loop level Feynman diagrams for b → sµ + µ − transitions. After this, the amplitudes are tensor decomposed in FeynCalc [134] and then, the resulting Passarino-Veltman functions are Taylor expanded in the external momenta up to second order. Finally, the functions are integrated with Package X [135]. In this way, with (3.12) where 7,8 defined in appendices C1 and C2 in [137]. Here, λ ii are the diagonal Yukawa couplings defined in SuperIso, G F is the Fermi constant and s W is the sine of the Weinberg angle. The Green function B H(0) for the box diagram contribution in C box 9,10 coming from the new lepton flavour violating (LFV) couplings is given by (3.17) Our computation shows two new terms absent in both the SuperIso manual and in [70], namely the mixed term in the C Z 9 expression where and a gauge dependent contribution to C box 9 coming from the box diagrams in figures 2c and 2d proportional to B . For all remaining terms, we obtained full agreement with [70] once the overall √ 2 factor in their Yukawa Lagrangian is taken into account compared to our Eq. (2.22). It is important to mention here that once the full quantum field theory matches with the effective theory at a scale µ W = O(m W ), the evolution of the WC C 7 (and C 7 ) from µ = µ W down to µ = µ b , where µ b is of the order of m b , is given at LO by [138] is given in Eq.(12.23) of [138] and references therein. The C 2 coefficient comes from four-quark operators generated by W boson exchange in the SM and contributes importantly when computing the branching ratio BR(B → X s γ). In the GTHDM, as shown in [70], an analogous contribution comes from charged Higgs exchange at tree level. In this way, following [139] with α S (m Z ) = 0.117, we use the following parametric expression at LO: where C 2 = C SM 2 + ∆C 2 for C SM 2 = 1 and with g 2 the weak coupling constant. Similarly, there will be a contribution to the C 9 (and C 9 ) WC coming from those four-quark operators given by [70] which can be added at LO to both the penguins and boxes contributions, obtaining

Summary of contributions
As already mentioned, in view of the flavour changing couplings in the GTHDM, there are two new contributions compared to the ones present in SuperIso. These contributions come from the box diagrams in figures 2a-2b and from the Z penguin in figure 1. The Γ L tb Γ L ts contribution is the largest and dominates the amplitude for most of the parameter space, with a strong dependence on tan β, m H ± , Y u 2,ct/tc , Y u 2,tt , Y l 2,µµ , Y l 2,µτ . There are also two subdominant contributions, the first one coming from the part proportional to Γ R tb Γ L ts in the Z penguin diagram in figure 1. When comparing its contribution relative to the Γ L tb Γ L ts term, we find regions of the parameter space in which it can make up to 10% of the total contribution (see figure 3 left). The second subdominant contribution is the already mentioned gauge dependent part of the boxes diagrams (figures 2c-2d) which is suppressed by the muon mass (see figure 3 right). Additionally we verified that when varying the mass of the charged Higgs from 500 GeV to 4000 GeV these ratios were essentially unaffected. In this way, we keep in our calculations the Γ R tb Γ L ts term from the Z penguin and neglect the gauge dependent part of the boxes diagrams. and C box,mix 9 refers to the first and second terms in Eq. (3.12) respectively.

b → c ν semileptonic transitions
As a consequence of the new interactions between the fermions and the charged Higgs, semileptonic tree level flavour changing transitions appear in the GTHDM (figure 4) which have been extensively studied in the literature [15,16,127,[140][141][142]. Therefore we include tree-level calculations of the Wilson coefficients related to these in our analysis. The effective Hamiltonian responsible for the b → c ν transitions for the semileptonic decays of B-mesons, including the SM and tree level GTHDM contributions can be written in terms of scalar operators in the form where C cb SM = 4G F V cb / √ 2 and the operators are given by Given that the flavour of the neutrino in the final state can not be discerned by experiments, one has to add (incoherently) to the SM the NP contributions associated with the LFV couplings ξ l ij . As the existing constraints will apply separately to the scalar and the pseudoscalar couplings, it is convenient to define where in our analysis we evaluate the WCs C cb R and C cb L at tree-level, with the expressions, (3.28)

Observables
In this section we present the observables to be included in the fit. We divide them in four sets: The first one for FCNCs in b → s transitions and B meson rare decay observables, both of them affected by the new WC contributions. The second set is associated with FCCCs observables that arise from semileptonic b → c ν decays and the mass difference ∆M s from B s − B s oscillations. Various leptonic decays of mesons form the third set. Finally, the fourth set contains leptonic observables associated with τ and µ decays, among them the anomalous magnetic moment of the muon in particular.

FCNCs and B rare decays
Lepton flavour universality in the SM means that all couplings between leptons and gauge bosons are the same (up to mass differences). This implies that any departure from this identity could be a clear sign of NP. The most interesting tests of LFU violation with FCNC are given by the ratios of b → sll transitions with Γ representing the decay width and K ( * ) are kaons. As per our choice of Yukawa textures in Eq. (2.27), here we only consider NP effects coming from the muon specific WCs, i.e., electronic WCs are SM-like. Aside from this R(K ( * ) ) ratios, hints for LFU violation are found in many branching fractions and angular observables related to B → K ( * ) µ + µ − decays as a function of the dimuon mass squared q 2 . In this work we use the same observables as in [65], with the predicted values obtained with SuperIso and with likelihoods provided via HEPLike. In particular, among the observables included are the optimised angular observables P ( ) i which have been constructed in order to minimise the hadronic uncertainties emerging from form factor contributions to the B 0 → K * 0 µ + µ − decay at leading order [143]. In view of that, experimentally these observables are obtained by fitting q 2 -binned angular distributions and they are defined in the theory as CP-averages integrated in the q 2 bins: , (4.2) where the J i functions and the normalisation constant N bin are given in [65]. Additionally, they can be related to the form factor dependent observables S i [144] as where A FB is the forward-backward asymmetry of the dimuon system and F L is the fraction of longitudinal polarisation of the K * 0 meson. The most sensitive observable to scalar operators is the branching ratio BR(B s → µ + µ − ) which also depends on the muon specific C 10 and C 10 WCs [65]: where f Bs is the decay constant and τ Bs is the mean lifetime. With respect to the inclusive B → X s γ decay, we use the full expression given in the works of [145][146][147][148][149][150] and implemented in SuperIso. The WCs C 7 and C 7 are constrained by this decay, given at the quark level by b → sγ, which at LO is We also take into account the rare decays B s → τ + τ − and B + → K + τ + τ − as well as the LFV processes B s → µ ± τ ∓ , B + → K + µ ± τ ∓ and b → sνν with theoretical expressions given in [70]. A list of the included FCNC observables 4 can be found in Table 1.
Observable Experiment Table 1: Experimental measurements of FCNCs observables and bounds for rare B decays considered in our study. The R νν K ( * ) parameters are related to b → sνν transitions as introduced in Eq.(4.6) in [70]. We also include all the angular distributions and branching fractions of B 0 → K * 0 µ + µ − decays, the branching fractions of both B s → φµ + µ − and B + → K + µ + µ − with measurements provided by the HEPLikeData repository [157].

FCCCs observables
The most relevant FCCC observables are the ratios of semileptonic B meson decays to τ and light leptons, that is where D ( * ) are charmed mesons and l is either an electron (e) or a muon (µ). As of the time of writing, the world average for the experimental measurement of the ratios R(D ( * ) ) sits at a 3.1σ deviation from the SM prediction [62]. The GTHDM contributions to R(D) and R(D * ) from the effective Hamiltonian in Eq. (3.25) can be written as, . (4.10) 4 New measurements of BR(Bs → µ + µ − ) have been performed recently by LHCb [151,152], as well as a combination with previous results [64], giving a combined measured value of 2.85 +0.34 −0.31 . Nevertheless, we do not expect significant deviations from our results with this new measurement.
In addition to R(D) and R(D * ), a third ratio has been measured by the Belle collaboration [158], the ratio R e/µ = BR(B → Deν)/BR(B → Dµν) which is considered to be the stringent test of LFU in B decays. It can be expressed in the GTHDM as where we have obtained the NP leptonic contributions by integrating the heavy quark effective theory (HQET) amplitudes of the scalar type operators from [159,160]. The B c meson lifetime has contributions from the SM, given by τ SM ps [161], and the GTHDM, which can be written as where the -1 term accounts for the subtraction of the SM contribution. By using the lifetime of the B c meson as the constraining observable, we can compare it to the current experimental measurement of τ Bc = 0.510 ± 0.009(ps) [155], instead of using the theoretical limits on the branching ratio BR(B c → τν), which are reported to be either 10% [162] and 30% [142] 5 . Another related measurement, B + c → J/ψτ + ν τ , has been reported by LHCb [166] and also hints to disagreement with the SM. However the errors are too large at present to reach a definitive conclusion, with R(J/ψ) = 0.71 ± 0.17 ± 0.18. In addition it has been claimed that the hadronic uncertainties are not at the same level as for the observables related to B → D * transitions [159], so we do not include it in our fit.
In contrast, a recent measurement of the longitudinal polarization fraction of the D * meson, defined as has been recently announced by the Belle collaboration [167], deviating from the SM prediction F SM L (D * ) = 0.457 ± 0.010 [168] by 1.6σ. The B → D * τ ν differential decay width into longitudinally-polarized (λ D * = 0) D * mesons is given (keeping NP from scalar contributions only) by where the helicity amplitudes are defined in appendix B of [159]. In addition, we also include the normalised distributions dΓ(B → Dτ ν)/(Γdq 2 ) and dΓ(B → D τ ν)/(Γdq 2 ), as measured by the BaBar collaboration [169]. Lastly, the mass difference ∆M s of B s − B s oscillations is included in our study and (for m A = m H ) is given by [71] ∆M GTHDM Bs = 1.022 [170,171], and the U running matrix being defined in [71]. A summary of all FCCC observables included in this study is provided in Table 2.

Leptonic decays of mesons
Beyond those described in Sections 4.1 and 4.2, there are additional leptonic decays included in this study. The total decay width at LO for the process M → lν in the GTHDM is computed as [126,132,172] where i, j are the valence quarks of the meson M , f M is its decay constant and ∆ ll ij is the NP correction given by where the relations depend on the Yukawa textures. The list of fully leptonic decays of mesons included in this analysis, for various mesons M , can be seen in Table 3.

Leptonic observables
There are a number of leptonic processes that are forbidden or suppressed in the SM but can occur in the GTHDM. These include modifications to the form factors for γ, Z and other interactions, which lead to contributions to the anomalous magnetic moment of the muon, (g − 2) µ , and LFV decays such as τ → µγ, τ → 3µ and h → τ µ. In the SM, the contributions to these LFV observables are suppressed by the GIM mechanism, giving a very low experimental background, but in the GTHDM LFV is allowed at oneand two-loop level through the couplings ξ l ij in Eqs. (2.24-2.26,B.3). 6 A second Higgs doublet has been examined as a way to explain the muon g −2 anomaly. In the Type-X [18][19][20][21][22][23][24][25][26][27][28][29][30][31] and Flavour-Aligned [32][33][34][35][36][37] versions of the THDM the contributions from two-loop diagrams are dominant in most of the parameter space thanks to mechanisms also available in the GTHDM. Additionally, with LFV, the one-loop diagrams can receive a chirality flip enhancement from including the tau lepton in the diagram loop, as was investigated by [38][39][40][41][43][44][45][46], however they only examined muon g − 2 contributions at the one-loop level.
Due to the similarity of the diagrams between → γ and muon g − 2 (which is effectively µ → µγ, see figure 5), these two observables share nomenclature and contributions. For both muon g −2 and τ → µγ we can break the contributions into the same three groups: one-loop, A ijL,R ; two-loop fermionic, A (2,f ) ijL,R ; and two-loop bosonic, A (2,b) ijL,R , contributions, so that the observables can be written as τ µL,R and α EM is the fine structure constant. All form factors A (l) ijL,R have been appropriately renormalised by combining with the relevant counterterms, and are all calculated using masses and couplings that have been extracted from data at tree-level. Additionally, for the contributions to muon g − 2 we must subtract 6 Note that in this study we will focus solely on the decays involving τ and µ leptons due to our choice of including only second and third generations in the ξ l ij matrix from Eq. (2.27).  off the SM contributions from the SM Higgs boson to obtain a purely BSM contribution to muon g − 2.
The entire one loop contribution for muon g − 2 and → γ can be found by summing over the neutral scalars φ and lepton generations: where the functions A  figure 5. To obtain the BSM contributions to muon g − 2, we must also subtract off the contribution from the SM Higgs boson to obtain a truly-BSM one-loop contribution. At the two-loop level we consider the Barr-Zee diagrams, shown in figures 6 and 7. Just as for the one-loop contributions before, we can subdivide each of these contributions The total fermionic two-loop contribution to muon g − 2 is given by [34] A  Note that these contributions do not include 2-loop diagrams with an internal Z boson leg, as in [32].
In the case of the τ → µγ decay, the contributions from the fermionic and bosonic Barr-Zee two loop diagrams, A  7 We do not consider two-loop bosonic diagrams that are not Barr-Zee diagrams, since their maximum contributions to muon g − 2 are relatively small [35], whereas Barr-Zee contributions have been proved to be dominant for some regions of the parameter space [176]. Additionally, two-loop diagrams involving neutral bosons where both legs are Higgs bosons are suppressed by a factor m 4 µ , while diagrams with both legs being either γ or Z are SM contributions, so we do not consider either, only those with both a φ and a γ or Z boson leg. Similarly for diagrams involving charged legs of H ± ,W ± , we only consider a H ± and W ± boson paired together, as a pair of H ± legs lead to diagrams with suppressed contributions [32].
The contributions to τ → 3µ decay can be divided up into 3 separate groups, the tree-level, dipole, and the contact contributions. The contributions from tree-level decay are computed in [127]. We have found that the dipole contributions, which involve the penguin-photon diagrams of the form of τ → µγ decays, are quite sizable compared to those at tree-level and cannot be ignored. Namely, they are given by [177]: Similarly, the contact terms involving effective four-fermion interactions [178] could be at first comparable to the dipole contributions. The contact contributions are given by where the coefficients g 2 and g 4 are given in appendix B. Another observable that we include is the lepton violating h → τ µ decay. This is given at tree level by 8 with the total decay width of h given by Γ h = 3.2 MeV [175]. Lastly, besides g − 2 and LFV observables, experiments have also provided constraints for the LFU ratio in τ decays. This ratio is commonly known as (g µ /g e ) 2 and is given as [126,172] where f (x) = 1 − 8x + 8x 3 − x 4 − 12x 2 log x and R ij is the BSM scalar contribution, given in the GTHDM as All of the experimental measurements and upper bounds for leptonic observables are shown in Table 4. 8 We computed the contributions coming from one-loop diagrams with two charged Higgses in the loop and found them to be 7 orders of magnitude suppressed compared to the tree level. Diagrams involving a pair of heavy neutral Higgses are possible as well but even more suppressed. The GTHDM only takes into account the tree level, which relies on being close to the alignment limit but not exactly, otherwise this tree level contribution would be zero.

Results
Our main goal is to study the impact of these observables on the GTHDM parameter space and, in particular, infer the goodness-of-fit of the model in light of these anomalies. Given the plethora of observables defined in the previous section and the large multidimensional parameter space, it is very important to combine them in a statistically rigorous manner in a global fit. This avoids serious shortcomings from more naive approaches like simply overlaying constraints from confidence intervals [181].
To visualize the results we will project the high dimensional parameter space onto twodimensional planes. To this end, the central quantity of interest is the profile likelihood, which is, for fixed parameters of interest θ 1 and θ 2 , the maximum value of the log-likelihood function that can be obtained when maximizing over the remaining parameters η. All profile likelihood figures in this study are created with pippi [182].
As mentioned earlier, we use here the GAMBIT framework for our study. The theoretical predictions of the model and the experimental likelihoods are either implemented natively in GAMBIT or from external tools interfaced with GAMBIT. In particular, the likelihoods related to b → sµ + µ − transitions are obtained from HEPLike, which retrieves experimental results and their correlated uncertainties from the HEPLikeData repository. To efficiently explore the parameter space, we employ the differential evolution sampler Diver, which is a self-adaptive sampler. We choose a population size of NP = 20000 and a convergence threshold of convthresh = 10 −6 . The data we present in this work comes from scans that took between 6 and 8 hours of running time on the Australian supercomputer GADI with cores varying between 1400 and 2000.

Parameter space
We perform the parameter scans in the physical basis, i.e., where the tree-level masses of the heavy Higgses, m H , m A and m H ± are taken as input. The remaining model parameters are tan β, m 12 and the Yukawa couplings Y f 2,ij as in Eq. (2.23). In order to avoid collider constraints, we work in the alignment limit choosing s β−α close to 1, and we select a conservative lower limit on the masses of the heavy Higgses m H,A,H ± ≥ 500 GeV 9 . We also fix m A = m H in our study, motivated by the requirement to satisfy the oblique parameter constraints which favour small mass splittings and in order to simplify the sampling of the parameter space. So as to choose reasonable priors for the Yukawa couplings, we take into account various constraints on them (or equivalently on ξ f ij ) from previous studies. The tighter theoretical constraints come from perturbativity which requires ξ f ij ≤ √ 4π ∼ 3.5.
On the phenomenological side, the studies in [43,131] have found values as large as ξ u tt ∼ 0.1 and ξ u tc ∼ 0.32 for masses of the heavy Higgses of order 500 GeV. With respect to the ξ u cc coupling, it has been shown in [132] that O(1) values are possible within the charged anomalies, and similar values were considered in [70] in the context of the neutral anomalies, not only for ξ u cc but for all the Yukawa matrix elements. As for the new leptonic couplings, the results in [15,70,176] indicate they should be O(1) or less in order to fit the charged anomalies. Lastly, the extra down Yukawa couplings ξ d ij are in general expected to be O(0.1) [40,184] and in particular ξ d sb is expected to be strongly constrained by B s − B s mixing. 9 From preliminary results we found that low Higgs masses are disfavoured by the contribution of various constraints and thus we do not attempt to include precise constraints on the masses from BSM Higgs searches (see e.g. [183] for a discussion of the limits on the charged Higgs mass). We leave a detailed collider study to future work.
The results of our scans show two degenerate regions of solutions according to the sign of Y u 2,tc . We indeed verified that these regions are degenerate and the final results are unaffected by this choice, hence we select Y u 2,tc ∈ [−2, 0] for the phenomenological analysis from now on. Namely, this degeneracy is a result of the dependency of various observables on products like Y u 2,tc Y f 2,ij where Y f 2,ij also flips its sign. 10 figure 8 we observe that the preferred values for the charged Higgs mass are of order 3 TeV with tan β ≈ 1. We find that the combined contribution of FCNC likelihoods fits better the data for this particular mass range. Similarly, although values of tan β up to 50 are possible in the GTHDM when using theoretical constraints alone, we identified that once we take into account all flavour constraints, there is a clear preference for low values, close to tan β ≈ 1, in agreement with [14,185,186]. This preference can be understood as follows. The box contributions in figures 2a-2b depend on the Green function B H(0) in Eq. (3.17), which for values of the charged Higgs mass m H ± < 2 TeV or m H ± > 4 TeV significantly over-or undershoot, respectively, the observed value of ∆C 9 ≈ −1 (see below). Furthermore, the measurement of the B c lifetime and the BaBar collaboration B → D ( ) τ ν distributions, both of which depend strongly on tan β and m H ± through the C cb R,L in Eq. (3.28), push both tan β and m H ± to values lower than 2 and greater than 2 TeV respectively. In addition to this, we have also noticed a strong penalty for large tan β values coming from the B s → µ + µ − decays, which is due to the strong tan β dependence on the C 10 and (pseudo) scalar WCs. Lastly the preferred masses of the other heavy Higgses, m H and m A , are of the same order as m H ± as was expected because of the oblique parameter constraints. The best fit values for some relevant scan parameters can be found in table 5. 10 We first found those two regions of solutions via an auxiliary scanning method based on the quadratic approximation to χ 2 as a function of the WCs (see appendix C).

Parameter
Best fit

Wilson coefficient
Best fit 0.002 ± 0.01 Re(∆C 7 ) 0.01 ± 0.01 Re(∆C 8 ) 0.002 ± 0.015 Re(∆C 8 ) 0.01 ± 0.01 Re(∆C 9 ) −0.89 ± 0.15 Re(∆C 10 ) −0.19 ± 0.14 Table 5: Best fit values for the scan parameters (left) and WCs for b → sµ + µ − transitions (right). We show only Re(∆C Q 1 ) given that at tree level and in the alignment limit Re(∆C Q 1 ) = Re(∆C Q 2 ) and m s /m b Re(∆C Q 1 ) = Re(∆C Q 1 ) = −Re(∆C Q 2 ). The uncertainties on the WCs were computed with GAMBIT assuming a symmetric Gaussian distribution from the resulting one-dimensional profile likelihoods. We do not display the Re(∆C 9,10 ) WCs either which we find to be suppressed by a factor of m b /m t compared to their non prime counterparts.

Neutral and charged anomalies
In table 5 we show the best fit values for the parameters from the scans (left) and the muon specific WCs evaluated at the best fit point (right), where in particular, ∆C 9 is consistent with the value obtained by model independent fits at the 1σ level. In this sense, the neutral anomalies can indeed be explained in the GTHDM as shown in figure 9. Furthermore, coming from the quadratic dependence in the branching ratio BR(B s → µ + µ − ), we can see two regions of solutions for the scalar WC ∆C Q 1 , one of them containing the SM prediction within 2σ. In addition, we ran a complementary scan invalidating points for |∆C Q 1 | > 0.1 and found that the corresponding region of solutions gives an equally good fit to the data, i.e., the preference over the second region of solutions is completely arbitrary. In order to better understand the contribution of the GTHDM to the various rates and angular observables, we display various plots comparing both the SM and the GTHDM predictions along the experimental data. For the angular observables P 1 and P 5 defined in Eqs.(4.4) and (4.5), we show in figure 10 their predictions compared to the CMS 2017 [187], ATLAS 2018 [188] and LHCb 2020 [189] data. For P 1 (figure 10 left) the GTHDM distribution is rather indistinguishable from the SM one, except in the [1, 2] GeV 2 bin close to the photon pole and sensitive to C ( ) 7 . The situation is different for P 5 (figure 10 right) in which the GTHDM prediction fits the LHCb 2020 data better, particularly in the C ( ) 7 -C ( ) 9 interference region (1 < q 2 < 6 GeV 2 ). We also provide in figure 11 predictions for the angular observables in the S i basis using the same LHCb 2020 measurements and also the ATLAS 2018 [188] data. We can see that the GTHDM fits better the LHCb data [189] in the large recoil region than the SM by 2σ. We also note that neither the SM or the GTHDM can explain the central values (with larger uncertainties) from the ATLAS 2018 data.
As for the measured branching ratios of B 0 → K 0 * µ + µ − and B + → K + µ + µ − , in figure  12 we show the SM and GTHDM predictions using the LHCb results [190][191][192], where we can see again how the GTHDM fits better the data compared to the SM, specially in the region above the open charm threshold, sensitive to both C ( ) 9 and C ( ) 10 . In contrast, the performance of the model is worse than the SM (figure 13 left) in the low recoil region of the differential branching ratio dBR dq 2 (Λ b → Λµ + µ − ) when comparing to the LHCb 2015 [193] data. As pointed out in [65], the decays of the Λ b baryon, such as Λ b → Λµ + µ − have much larger uncertainties than those of the corresponding meson decays. However, once more experimental data is available, recent [194] and future developments of lattice  calculations would eventually make this decay providing similar constraints as other b → sµ + µ − transitions. Finally, the results for the dBR dq 2 (B s → φµ + µ − ) distribution are shown in figure 13 right. The large recoil region of the experimental data deviates from both the SM and GTHDM predictions by approximately 3σ, and for the low recoil bin the GTHDM performs slightly better than the SM by approximately 1σ.  Figure 12: Left: Differential branching ratio for dBR dq 2 (B 0 → K * 0 µ + µ − ) with the LHCb 2016 data [190]. Right: dBR dq 2 (B + → K + µ + µ − ) compared to the LHCb 2012 and 2014 measurements [191,192].  [193] data. Right: dBR dq 2 (B s → φµ + µ − ) compared to the LHCb 2015 and 2021 data [195,196].
Last but not least important observables related to the b → sµ + µ − transitions are the ratios R(K ( * ) ). Despite being only three bins in total [151,153,197], these measurements have been intensively studied as they provide evidence for LFU violation. We include in our fit the latest LHCb collaboration data for the R(K * ) and R(K) ratios from 2021 [151] and 2017 [153] respectively and obtain the plots in figure 14, where we compare also to the Belle 2019 experiment data [198,199]. The effect from the fit on the R(K ( * ) ) ratios is significant, explaining the LHCb 2021 measurement of R(K) at the 1σ level.
The next interesting results are related with the charged anomalies, in particular we find that the R(D ( * ) ) ratio can (can not) be explained at the 1σ level with the GTHDM, a result in agreement with the phenomenological analysis of [132]. We furthermore corroborate that  Figure 14: R(K ( * ) ) theoretical ratios compared to both the LHCb [151,153] and Belle data [198,199]. the constraint coming from the B c lifetime makes it very difficult to fit R(D * ) and R(D) simultaneously. In figure 15 we show the preferred values by the profile likelihood. We see a slightly better performance of the GTHDM compared to the SM with respect to the HFLAV average. Regarding the dΓ(B → D ( ) τ ν)/(Γdq 2 ) distributions measured by BaBar [169], we find that the GTHDM prediction is indistinguishable from the SM, in agreement with [16]. We find furthermore that the longitudinal polarisation F L (D * ) is strongly correlated with R(D * ) and the model is not able to explain the Belle measurement, giving a best fit value of 0.458 ± 0.006.

Anomalous (g − 2) µ
With regards to the anomalous magnetic moment of the muon, (g − 2) µ , we find that a simultaneous explanation using all the likelihoods defined before is not possible (solid red line in figure 16). However, when doing a fit to all other observables except the neutral anomalies, i.e., without using the HEPLike likelihoods, the model is able to explain the measured ∆a µ by Fermilab at the 1σ level (dashed gray line in figure 16). Furthermore, when evaluating the performance of the HEPLike likelihoods for the best fit value, we find a SM-like behavior with all NP WCs close to zero, except for those scalar WCs that enter in BR(B s → µ + µ − ).

Projections for future and planned experiments
Although a detailed collider analysis is beyond the scope of the present work, we have included as pure observables the branching ratio for t → c h and h → b s 11 at tree level. These tree level branching ratios in the GTHDM are suppressed as c 2 βα |ξ u(d) tc(bs) | 2 , respectively, so that in the alignment limit they will be exactly zero. In order to study the effects of this fined tuned suppression, we have ran a second scan with s βα ∈ [0.9999, 1] and we found the branching ratio of t → c h decays are of order 10 −11 − 10 −7 , which although are outside future searches sensitivities, they are larger than the SM loop prediction (∼ 10 −15 ) and well below the current experimental upper bound obtained by the ATLAS collaboration [200] BR(t → c h) < 1.1 · 10 −3 .
Concerning the BR(h → b s) observable, it was shown in [71] that it is related to tree level B s − B s oscillations which are not only proportional to c 2 βα but also to pseudoscalar contributions independent of the scalar CP-even mixing. Hence, in figure 17 we see that h → b s is not as constrained as t → c h with values ranging from 10 −7 up to 10 −3 at the 1σ level, which may be within range of the ILC [201].
Regarding LFV searches, we show in figure 18 the profile likelihood for the τ → 3µ and τ → µγ branching ratios. We see that the best fit value for the τ → 3µ decay is well 11 We are not aware of current bounds for the h → b s branching ratio so we did not define an associated likelihood function for it.  [71].
within the projected sensitivity in the Belle II experiment [202] with a discovery potential for BR(τ → 3µ) ∼ 10 −9 . Regarding the τ → µγ decay, we find that with the projected future sensitivity, the GTHDM prediction could be confirmed with values for the branching ratio varying from 10 −9 up to 10 −8 . As mentioned earlier, the τ → 3µ decay receives contributions in the GTHDM from all tree, dipole and contact terms, in such a way that a possible detection in the τ → µγ channel will not necessarily imply a strong constraint for BR(τ → 3µ).  With respect to h → τ µ, with the model best fit point values, we computed the branching ratio BR(h → τ µ) obtaining values from 10 −2 down to 10 −6 which are within the future sensitivity at the HL-LHC, reaching the 0.05% limit [203].
Finally, as for the B s → τ + τ − decay, we find values of at most BR(B s → τ + τ − ) ∼ O(10 −6 ) with our best fit point, which is one order of magnitude higher than the SM prediction, but out of reach of the future sensitivity in both the HL-LHC and the Belle-II experiments with limits at O(10 −4 ) [202,204]. Regarding the branching ratio BR(B + → K + τ + τ − ), as in the B s → τ + τ − decay, the predicted branching ratio BR(B + → K + τ + τ − ) is of order 10 −7 − 10 −6 , out of reach for Belle-II projections at 2 × 10 −5 .

Conclusions and Outlook
We presented a frequentist inspired likelihood analysis for the GTHDM including the charged anomalies, b → sµ + µ − transitions and the anomalous magnetic moment of the muon along with other flavour observables. The analysis was carried out using the open source global fitting framework GAMBIT. We computed the GTHDM WCs and validated them obtaining full agreement with the one loop calculations reported in the literature after the different notation factors were taken into account. As expected, we found that the GTHDM can explain the neutral anomalies at the 1σ level. Additionally, we also confirmed that the model is able to fit the current experimental values of the R(D) ratio at the 1σ level, but it can not accommodate the D * charmed meson observables R(D * ) and F L (D * ). Furthermore, we inspected the fitted values for the angular observables in b → sµ + µ − transitions, obtaining in general a better performance with the GTHDM in comparison to the SM.
Then, based on the obtained best fit values of the model parameters and their 1σ and 2σ C.L. regions, we made predictions impacting directly in the future collider observables BR(t → ch), BR(h → bs), BR(h → τ µ), BR(B s → τ + τ − ), BR(B + → K + τ + τ − ) and the flavour violating decays of the τ lepton, BR(τ → 3µ) and BR(τ → µγ). We find that the model predicts values of BR(t → ch), BR(B s → τ + τ − ) and BR(B + → K + τ + τ − ) that are out of reach of future experiments, but its predictions for BR(h → bs) and BR(h → τ µ) are within the future sensitivity of the HL-LHC or the ILC. We also find that the predictions for the τ → 3µ and τ → µγ decays are well within the projected limits of the Belle II experiment. In summary, the next generation of particle colliders will have the sensitivity to probe, discover or exclude large parts of the parameter space of the GTHDM, and thus it serves as a further motivation for the development of higher energy and higher intensity particle colliders.
We can envision many avenues of future investigation using the tools and techniques developed for this work. The complete parameter space of the GTHDM is enormous, and thus for this study we have only focused on a subset of CP-conserving interactions between second and third generation fermions. The inclusion of the first generation in the Yukawa textures would introduce additional interactions and decay channels, possibly improving the fit to various of the flavour anomalies, while at the same time introducing new relevant constraints, such as rare kaon decays, e.g. from the NA62 experiment, and LFV muon decays, e.g. from the Mu2e experiment. CP-violation in kaon and B-meson decays would also become important constraints in case of complex Yukawa textures. Modifications of the GTHDM may also lead to improved fits to some flavour observables, for instance it has been shown that with the addition of right-handed neutrinos the model can better accommodate the neutral anomalies. Lastly, in this study we have not included detailed collider constraints from, e.g., searches for heavy Higgs bosons at colliders. Such a detailed study is a clear follow up from this work and it will showcase the complementarity of flavour and collider searches to constrain models of new physics that tools such as GAMBIT can explore.
Finally, in view of the latest experimental measurement made by the Fermilab Muon g − 2 Collaboration, we performed a simultaneous fit to ∆a µ constrained by the charged anomalies finding solutions at the 1σ level. Once the neutral anomalies are included, however, a simultaneous explanation is unfeasible. A detailed study looking for a simultaneous explanation of both g − 2 and the neutral anomalies in the GTHDM will be presented in a follow-up work. where

B Loop Functions and Vertex Couplings
The one-loop contributions can be separated into fermion-fermion-scalar (FFS) and scalarscalar-fermion (SSF) diagrams shown in the diagrams shown in figure 5. As seen in these diagrams, we can have any one of the SM charged leptons or neutrinos paired with any neutral Higgs boson φ = h, H, A or the charged Higgs boson H ± respectively. The contributions from each of these diagrams with a scalar-fermion pair from a lepton of generation a to a lepton of generation b are shown below: . Additionally, to get the BSM contributions to muon g − 2, one must subtract of the SM contribution. This contribution is obtained from Eq. (B.1), A (F F S) µµL,R (h SM , µ), by using a mass of m h SM = 125.09 GeV, and a Γ f h SM coupling in Eq. (2.24) with c βα = 0 to obtain a SM-like coupling. For muon g − 2, the dominant BSM contributions at the one-loop level come from the chirality flipping term involves an internal τ lepton, m µ m τ /(48π 2 m 2 φ )Γ φτ µ Γ l * φµτ F m 2 τ /m 2 φ , with an enhancement of m 2 τ /m 2 µ . The coupling ξ τ µ should be nonzero to get the chirality flip enhancement from the internal τ lepton in figure 5.
The GTHDM contributions to muon g − 2 from Barr-Zee diagrams with a fermionic loop are given by [34] where c 2 W = cos 2 θ W , g f v = T 3f − 2Q f s 2 W , and T 3f denotes the isospin of the loop fermion (T 3f = (1/2, −1/2, −1/2) for (u, d, l)). The contributions A    where the masses have been ordered so that m 1 < m 2 < m 3 . The contributions to muon g − 2 from Barr-Zee diagrams with a bosonic loop are given by [32] W ∓ terms: (B.23) The loop functions used for the two-loop muon g − 2 bosonic Barr-Zee contributions come from [32] and are defined as one dimensional integrals: .

(B.26)
Similarly, the Barr-Zee fermionic and bosonic contributions to l → l γ are given by These do include the Z boson contributions as per [176]. The coupling g φW ± W ∓ is defined in Eq. (B.21), and t 2 W = tan 2 θ W .

C Auxiliary scanning method
The two regions of solutions for Y u 2,tc were expected already when applying the quadratic approximation to the χ 2 function defined in [206] for a fit to the b → sµ + µ − observables solely. Explicitly, the likelihood function is approximated by where C = {∆C 7 , ∆C 9 , ∆C 10 , ∆C 7 , ∆C 9 , ∆C 10 } are the WCs used as parameters to be fitted, and Cov −1 is the covariance matrix or Hessian obtained using the minuit and flavio packages, which encodes a fit using the likelihoods from [65] (excepting the associated likelihoods for the Belle experiment measurements not available in flavio). After obtaining the Hessian, a random generator in Mathematica is requested to find points inside the ellipsoid defined by ∆χ 2 ≤ σ 2d (1) and ∆χ 2 ≤ σ 2d (2) for 2 degrees of freedom and boundaries defined by the values of the parameter space in Eq. (5.2). With this auxiliary method, we were able to help the Diver sampler to scan over different regions of the parameter space.