Cornering the Two Higgs Doublet Model Type II

We perform a comprehensive study of the allowed parameter space of the Two Higgs Doublet Model of Type II (2HDM-II). Using the theoretical framework flavio we combine the most recent flavour, collider and electroweak precision observables with theoretical constraints to obtain bounds on the mass spectrum of the theory. In particular we find that the 2HDM-II fits the data slightly better than the Standard Model (SM) with best fit values of the heavy Higgs masses around 2 TeV and a value of $\tan \beta \approx 4$. Moreover, we conclude that the wrong-sign limit is disfavoured by Higgs signal strengths and excluded by the global fit by more than five standard deviations and potential deviations from the alignment limit can only be tiny. Finally we test the consequences of our study on electroweak baryogenesis via the program package BSMPT and we find that the allowed parameter space strongly discourages a strong first order phase transition within the 2HDM-II.


Introduction
Two Higgs Doublet Models (2HDM) [1] (see Ref. [2] for a textbook and Ref. [3] for a more recent review) are very well-studied extensions of the Standard Model (SM) and probably one of the simplest possibilities for beyond the SM (BSM) physics. Besides the SM particles, four new scalars arise in the 2HDM: two charged ones (H + and H − ), a pseudoscalar (A 0 , CP-odd) and a scalar (H 0 , CP-even; when denoting the SM Higgs with h 0 ). Such an extension might solve some of the conceptual problems that arise within the SM, see e.g. Ref. [4]. In particular, 2HDMs could provide the missing amount of CP violation to explain the matter-antimatter asymmetry in the Universe and they could also lead to a first order phase transition in the early Universe, as required by the Sakharov criteria [5]. 2HDMs further appear as subsets of more elaborate BSM models, such as the minimal supersymmetric SM (MSSM). Depending on the structure of the Yukawa couplings one distinguishes different types of 2HDMs -in this work we will consider the 2HDM-II.
Since we have no direct experimental evidence so far for an extended Higgs sector, we investigate indirect constraints on the 2HDM-II by comparing the most recent measurements with high precision calculations. We perform a comprehensive study of more than 250 observables, where the new Higgs particles could appear as virtual corrections and modify the SM prediction. Our study extends previous works like Refs. [6,7]. 1 In 2009, Ref. [6] studied bounds from quark flavour observables on the 2HDM-II within the framework of CKMfitter. In particular, it was found that the 2HDM does not perform better in the fit than the SM and that the new charged scalar has to be heavier than 316 GeV. We will considerably extend this analysis by using updated measurements, theory predictions and by including additional flavour observables. In 2017, the authors of Ref. [7] used HEPfit to constrain the 2HDM-II with Higgs data, electroweak precision observables and some flavour observables (the radiative decay b → sγ and the mass difference of neutral B mesons ∆m q ), as well as theoretical constraints like perturbativity. The main findings are lower mass bounds such as m H + > 740 GeV and a lower limit on the mass splittings of the heavy Higgses of the order of 100 GeV. Higgs signal strengths allow only small deviations from the alignment limit, |β − α − π/2| < 0.055 -this bound is further reduced to |β − α − π/2| < 0.02 within the overall fit (driven by theory constraints). Finally, the authors of Ref. [7] found that the parameter tan β lies in the preferred range [0.93; 5.0]. We will extend this analysis by updating experimental values and considering a large number of flavour observables. Moreover we confirm the observation made in Ref. [7], that the wrong-sign limit is excluded by data.
Throughout this work, we use the Python program package flavio [11] for flavour physics for our predictions of observables in the 2HDM-II. New physics models are not explicitly included in flavio. Instead, they impact observables through their contributions to the Wilson coefficients of dimension-6 operators in the relevant effective theory: the Weak Effective Theory (WET) below the electroweak scale, or the Standard Model Effective Field Theory (SMEFT) above. In Section 6, the relevant 2HDM-II contributions are included in the coefficients of the WET operators with five active flavours (there are only three active flavours in the WET-3 basis used for the anomalous magnetic moment of the muon a µ ). For a full list of operators in the flavio basis for each effective theory, see the documentation of the WCxf package [12].
Once 2HDM-II contributions to the relevant coefficients have been included, we use flavio to construct fits in 2HDM parameter spaces, where we show in colour the allowed regions of our parameter space within the confidence levels stated in the captions. The package flavio expresses experimental measurements and input parameters as various one-and two-dimensional probability distribution functions. The errors of theoretical predictions are constructed by calculating the prediction for each of a number (we use N = 10 4 ) of randomly-selected values of the input parameters within their probability distributions, and then computing the standard deviation of these values. We work with the "fast likelihood" method in flavio for constructing likelihood functions of the form L = e −χ 2 ( ξ)/2 , where it is assumed that the set of fit parameters ξ contributing to the observables entering the likelihood function are taken at their central values. This method uses the combined experimental and theoretical covariance matrix of the observables in the fit. For further information on the "fast likelihood" method and the package's construction, we refer to flavio's documentation [11] or the flavio webpage 2 .
The program package flavio is used in Section 6 for the flavour observables, where we work in two-dimensional parameter spaces. We make use of flavio to construct likelihoods and combine measurements in Section 5 for the Higgs signal strengths, however for the calculations of the signal strengths themselves, we use analytical expressions written in Python natively. The constraints discussed in Section 3 and 4 stemming from theoretical considerations and the electroweak precision observables cannot be expressed in terms of coefficients of dimension-6 operators. Therefore in Section 3, the results are taken from a native Python code generating Monte Carlo scans of the full parameter space. The likelihood functions for the electroweak precision observables, taking into account the experimental correlations were translated into native Python code to be properly combined in Section 7 with all other observables and theory constraints using some customisations to the flavio package.
The paper is organised as follows. Our notation for the 2HDM-II is set up in Section 2. In Section 3 theoretical constraints like perturbativity, vacuum stability, and unitarity are studied, while Section 4 introduces electroweak precision observables in the form of the oblique parameters. Constraints stemming from the SM Higgs production and decay at the LHC are studied in Section 5. Most of our observables stem from the quark flavour sector, see Section 6. In Section 6.1 we investigate leptonic and semi-leptonic tree-level decays, where the charged Higgs is predominantly responsible for the new contributions. In this section we further study a potential pollution of the determination of CKM elements from semi-leptonic and leptonic decays by the 2HDM-II. B-mixing is considered in Section 6.2 and loop-induced b quark decays are studied in Section 6.3. Section 6.3.1 contains a brief description of the b → sγ transition, known to give a lower bound on the mass of the charged Higgs boson. The leptonic decays B q → µ + µ − get also sizable contributions from the new neutral Higgses, see Section 6.3.2. In Section 6.3.3 we study semi-leptonic b → s + − -transitions, where the so-called flavour anomalies have been observed in recent years. These anomalies can be best explained by vector-like new effects, while we consider here the effects of new scalar couplings. The results of our global fit, combining all constraints discussed so far are presented in Section 7. In Section 7.2 we comment on the lepton flavour observable, a µ , the anomalous magnetic moment of the muon, where recent measurements at Fermilab [13] have confirmed the older BNL value [14]. We study two scenarios based on using the SM prediction from the theory initiative [15] or a recent lattice evaluation [16]. Using the program package BSMPT [17,18] we investigate in Section 7.3 the question of whether our allowed parameter space can also lead to a first order phase transition in the early Universe. Finally we conclude in Section 8, while all our inputs are collected in Appendix A.
In the SM, a conjugate of the Higgs doublet is used to ensure both up-type and down-type quarks acquire mass. This can also be achieved by using two distinct complex scalar doublets with i = 1, 2. These acquire different vacuum expectation values (VEVs) v i after electroweak symmetry breaking (EWSB), that are related to the SM Higgs VEV v by v 2 1 +v 2 2 = v 2 [1][2][3]. This construction has 8 degrees of freedom, 3 of which become the longitudinal polarisations of the massive W ± and Z bosons, leaving 5 physical Higgs bosons. Two of these are charged, H ± , two are neutral scalars, H 0 and h 0 , and one is a neutral pseudoscalar, A 0 . Typically h 0 is taken to be the lighter of the two neutral scalars. The potential for a general 2HDM with a softly broken Z 2 symmetry is, in the lambda basis [3], (2. 2) The Z 2 symmetry is required to forbid tree-level flavor changing neutral currents (FCNC). Constraints on the λ i from a theoretical perspective are considered in Section 3. The parameters m 2 11 and m 2 22 are related to v 1 and v 2 as [19] In the remainder of this work we make use of the mass basis, through a transformation [20] 3 + v 2 λ 1 cos 2 α cos 2 β + λ 2 sin 2 α sin 2 β + λ 3 + λ 4 + λ 5 2 sin 2α sin 2β , The mass basis allows us to concentrate on the 8 physical parameters: the masses of the new particles, m H ± , m h 0 , m H 0 , m A 0 , the mixing angles α and β, the softly Z 2 breaking term m 12 , and the VEV v.
We identify h 0 as the observed SM-like boson with a mass of 125.1 ± 0.14 GeV [21] and take v = 246 GeV throughout this work. The angles α and β, describing the mixing of the neutral and charged Higgs fields, respectively, satisfy the following relations [20] tan (2.10) One may revert Eqs. (2.5) -(2.10) and express explicitly the λ i parameters in terms of the physical parameters, see e.g. Ref. [22]. We will show below that data strongly favours the small cos(β − α) limit, where the expressions for the λ i become particularly simple [22,23]: The Lagrangian for the Higgs and Yukawa sectors in the 2HDM is given by where L Y depends on the details of the quark-Higgs couplings in the 2HDM. There are four main types of 2HDM, based on which particles each of the doublets couple to. In this work, the 2HDM of Type II (2HDM-II) is investigated as it gives quark mass terms in the Lagrangian that are similar to the form in the SM. This in turn yields a similar CKM matrix through which all flavour-changing interactions can be described, with tree-level FCNCs forbidden by the imposition of the Z 2 symmetry (which is only softly broken by m 2 12 ). The 2HDM-II model has the doublet Φ 1 coupling to down-type quarks and leptons, while Φ 2 is coupling to up-type quarks, so are the Yukawa couplings, i, j = 1, 2, 3 run over the generations of the fermions, and the left-handed SU (2)-doublets of quarks and leptons are, respectively, -5 -After EWSB, rotating to the mass basis and focusing on the gauge and Yukawa couplings of h 0 gives Lagrangian terms of a very similar form to the SM, with factors κ i relating the coupling constants to those in the SM [24]: 18) where, in the 2HDM-II, Naturally, in the SM all κ i = 1, which is consistent with LHC data [25][26][27]. Setting cos (β − α) = 0 recovers this and thus matches the phenomenology of h 0 in the 2HDM with the SM Higgs. As such, this is known as the alignment limit. The dependence of these, and numerous other, couplings on cos (β − α) and tan β leads us to choose these as parameters of the model, in place of α and β.
The addition of new Higgs fields introduces new interactions, where the charged Higgs fields can replace W ± fields in flavour-changing charged interactions, as well as direct couplings of the new bosons to the weak bosons themselves. The part of the Yukawa Lagrangian describing interactions of H ± bosons with the fermions is given by [2,3] 4 where P L,R = (1 ∓ γ 5 )/2, m u and m d are the masses of the up-and down-type quarks, m ( = e, µ, τ ) is the charged-lepton mass, and V ud denotes the corresponding element of the quark-mixing Cabibbo-Kobayashi-Maskawa (CKM) matrix.
The new Higgs fields, discussed above, will affect many measured particle physics quantities including the electroweak precision parameters, Higgs signal strengths and flavour physics observables. These effects are investigated in Sections 4, 5, and 6 respectively and used to set bounds on the parameter space of the model. However, we start with the discussion of purely theoretical constraints.

Perturbativity
Within the lambda basis, perturbativity in the scalar sector can be simply expressed as [28,29] Clearly this is not a strict experimental constraint like the bounds studied below; it is more related to our ability to make meaningful predictions in the 2HDM -and the numerical value of the bound on the λ i contains a certain amount of arbitrariness. Thus we will also consider the less conservative bound of |λ i | ≤ 4 -inspired by the results in Refs. [30,31]. None of our conclusions will actually be changed by modifying the perturbativity bound from 4π to 4; to be able to study some more extreme scenarios in Section 7.3 we will, however, use the more conservative bound 4π. Writing the expressions for λ i in terms of the additional Higgs masses we find that for low masses, i.e. M 1000 GeV, no bounds (below the particle mass M ) are set on the mass differences of the new Higgs particles, while for high masses a form of degeneracy has to hold. In our analysis, we will consider heavy Higgs masses ranging between 10 2.5 ≈ 300 GeV and 10 5 GeV. Masses of the charged Higgs boson as low as 300 GeV are clearly ruled out by flavour observables, in particular by b → sγ, as discussed below. The experimental lower bound on the charged Higgs mass from direct searches is actually only m H + 160 GeV [32]. The upper bound has been chosen to be equal to the centre-of-mass energy of a future 100 TeV collider.
Looking at the couplings in Eq. (2.20), one can derive two further perturbativity constraints from the Yukawa sector where we again have chosen a very conservative range. For the range of log(tan β) we will consider the conservative lower bound tan β = 10 −1.5 = 0.03, and the upper bound tan β = 10 +2.5 = 300.

Vacuum Stability and Unitarity
Next we apply the conditions for a stable vacuum as set out in Ref. [33]: Demanding the vacuum to be the global minimum of the potential, we further find [34] We also consider the conditions from tree-level unitarity, see Refs. [29,35] 5 , alongside NLO unitarity and the condition that NLO corrections to partial wave amplitudes are suppressed relative to LO contributions, see Refs. [30,31].
To implement these constraints alongside the perturbativity conditions, we carry out Monte Carlo scans. Starting with the allowed sets of λ i (via Eq. (3.1)) we then use Eqs. (2.6) and (2.10) to get a relation among the parameters λ 1 ,..., λ 5 , β and the small parameter δ := β − α − π/2 6 , which is explicitly given in Appendix B. Expanding the expressions up to second order in δ, Eqs. (2.5)- (2.8) give the allowed ranges of Higgs masses depending on the value of β and δ. We show the effect of the theory constraints on the Higgs masses m H + , m H 0 , m A 0 , tan β and cos(β − α) as two-dimensional projections in Fig. 1. The starting values |λ i | < 4π are depicted in green, while points from |λ i | < 4 are shown in amber -the latter obviously leads to more constrained regions. Here we note that due to the imposition of NLO unitarity constraints, for starting values |λ i | < 4π, we find most successful points for |λ i | 2π. Moreover, we find that for charged Higgs masses above 1 TeV -as implied by the flavour observables -the choice of |λ i | < 4π or |λ i | < 4 has only a very minor impact. Figure 1: Bounds on the heavy Higgs masses m H 0 , m A 0 and m H + , tan β and cos(β − α) stemming from the theory constraints (perturbativity, vacuum stability and unitarity conditions). For the plots 10 8 points were generated and only values of |δ| < 0.5 were considered. The starting values |λ i | < 4π are given in green, while |λ i | < 4 is shown in amber.
The two upper plots and the left of the middle plots of Fig. 1 show that theory constraints require the new Higgses to be more or less degenerate in mass. The heavier the mass, the stronger this requirement becomes, which can also be nicely read off from Table 1. The right plot of the middle row in Fig. 1 shows that |δ| will be constrained to small values, if the new Higgs particles are heavier than about a TeV (such a bound on the charged Higgs mass will follow from flavour constraints discussed below). Again, the constraints on |δ| become stronger the heavier the new Higgs masses become. Quantitative values for the bounds on |δ| are given in Table 1. In the lower row of Fig. 1 we see that tan β is not really constrained by the theory constraints, apart from the restrictions given in Eq. (3.2).
In the alignment limit, the value of α is fixed by the value of tan β, and we find a similar requirement for a mass degeneracy of the new Higgs particles as in the general case. Moreover, we do not find any bound on tan β within the regions of interest.  Table 1: A numerical summary of the effects of imposing the theory constraints with starting values |λ i | < 4π. These intervals are found by carrying out a 10 8 points Monte Carlo scan with |δ| < 0.5, m H ± , m A 0 , m H 0 < 5000 GeV and −1.5 < log 10 (tan β) < 2.5. We have defined ∆ ij := m i − m j .

Electroweak Precision Observables
The oblique parameters first proposed in Ref. [37] are defined as Π XY (q 2 ) represents the self energies of the electroweak gauge bosons X and Y with 4-momentum q µ , the subscript i represents the W i field and the subscript Q represents the B field. Here we stress that in defining these oblique parameters a subtraction of the SM contributions is made; S = S NP − S SM and similarly for T and U , such that, by construction, in the SM S = T = U = 0. While most of the constraints we consider below depend on the particular type of the 2HDM, the oblique parameters are, at one loop order, unaffected by the Yukawa couplings and therefore universal for all 2HDMs. The expressions for S, T , and U in the 2HDM (derived from Ref. [38]) have been collected in Appendix C.
New physics generally has only a small effect on U ; we therefore follow the approach outlined in Ref. [21] and set U to zero. This effectively reduces the error on the experimental result for T , due to the correlation between T and U . The oblique parameters will be included in our global fit.

Higgs Signal Strengths
Signal strengths are a key measurement in determining if the observed Higgs boson is indeed that of the SM, or if new physics is involved. For a single production mode i with a cross section σ i , and decay channel f with a branching fraction B f , the signal strength µ f i is defined as as the cross section and branching fraction cannot be separately measured without further assumptions. Evidently if the SM is accurate, all these signal strengths should take a value of 1.
In practice, the small branching fractions and cross sections render some channels impossible to measure currently, but there are recent measurements available for 31 channels from the ATLAS and CMS collaborations at √ s = 13 TeV with an integrated luminosity of up to 139 fb −1 and 137 fb −1 respectively [25][26][27][39][40][41][42]. These include various combinations of the four main h 0 production channels at the LHC: gluon-gluon fusion (ggF), vector boson fusion (VBF), in association with a weak vector boson (Vh) or a top-antitop pair (tth), and the decay channels γγ, ZZ * , W W * , Zγ, τ + τ − , bb and cc. The recent observations of h 0 → µ + µ − [42,43] are also included here. In cases such as h 0 → µ + µ − where a combined signal strength across all channels is given, it is interpreted as being purely from gluon-gluon fusion as this mechanism dominates the Higgs production cross section at the LHC. In instances where a correlation matrix between measurements is given this is also included in the analysis. For full details on which channels have been used and their experimental values, see Table 4 in Appendix A. Where multiple measurements are available for a single channel we take an average, which significantly improves our sensitivity (see flavio's documentation for details [11]). The couplings of h 0 to the massive gauge bosons and fermions are multiplied by factors κ i compared to the SM Higgs, as set out in Section 2, with the relations for κ i given in Eq. (2.19). As the couplings are modified, the signal strength data can be used to place constraints on the two parameters on which the κ i depend: cos (β − α) and tan β. A scenario in which has been introduced as one of the ways a fourth sequential chiral generation of fermions may still exist 7 while remaining hidden from current LHC data [46]. By enforcing κ u = −κ d = 1 one finds Simply adding a fourth sequential, chiral, perturbative generation of fermions, without any additional new degrees of freedom is definitely ruled out, see e.g. Refs. [44,45]. The full criteria of the wrong sign limit are satisfied exactly only when this equation is followed at very large tan β, giving the set of couplings in Eq. (5.2). The contributions of the additional Higgs bosons to the loop-level processes have been neglected as they are negligible for the allowed masses found in Section 6.
The results of the analysis including all the Higgs signal strengths listed in Table 4 are shown in Fig. 2 (see also Table 2); we will later include the Higgs signal strengths in our global fit. This fit is performed using the analytic expressions for the Higgs production and decay channels in the 2HDM found in Ref. [2]. As outlined above, the coupling modifications, the κ i , are functions of tan β and cos(β − α), allowing us to constrain this plane. For some recent studies that include a similar analysis in the 2HDM, see, for instance, Refs. [7,23,24,47,48]. We observe that the alignment limit must be closely followed, which is unsurprising given that almost all of the experimental data included here matches the SM within 2σ, with the majority agreeing at 1σ. As a result, the 2HDM must closely replicate the SM behaviour for h 0 . At 2σ we find with the maximum allowed value of cos (β − α) occurring at tan β ≈ 1. Due to the dependence of κ d, on tan β, this maximum value falls significantly once away from tan β ≈ 1, to the point where, at large values of tan β, the alignment limit must be strictly followed and it is valid to set cos (β − α) = 0. At 3σ and above we find the wrong sign limit to be allowed, which substantially increases this upper bound, provided Eq. (5.3) is closely followed. : Contour plot of the coupling constant modifiers, in which we have fixed to its best fit value κ V = 1.036 and marked the alignment (wrong-sign) limit in black (red). The contours indicate the allowed space at 1, 2, 3, 4, 5σ from dark to lighter.
The wrong sign condition from enforcing κ u = −κ d = 1, Eq. (5.3), which is shown in Fig. 2 as a dashed red line, is here found to be excluded at 2.7σ by the signal strength data, an improvement over previous works such as Refs. [23,25,49], due in part to our inclusion of a larger data set than in such works, which improves our sensitivity.
To further probe the wrong-sign limit, we examine the possible values of the coupling modification factors, independently from Eq. (2.19). We first perform a fit in terms of the three κ i to find the best fit point, and then fix one of the κ i to its best fit value, allowing a 2D fit in terms of the remaining κ j,k =i to be performed. We show the results from such a fit in the (κ u − κ d, ) plane in Fig. 3, where we have fixed to the best fit value κ V = 1.036, and stress again that these results do not make use of Eq. (2.19). This fit does not exclude κ d = −1, in line with the current lack of sensitivity to the sign of κ d at the LHC. We then look to find solutions in terms of tan β and cos(β − α) that give values of κ u and κ d that lie within the 2σ allowed region close to κ d = −1 (bottom part of Fig. 3) by using Eq. (2.19). However, the only solutions we find require that | cos(β − α)| 0.96. Such extreme values of cos(β − α) are not consistent, in the 2HDM-II, with the fixed value of κ V = 1.036 from the best fit point that we use in the scan. Therefore this 2σ allowed region in the (κ u − κ d, ) plane cannot be attained in the 2HDM-II. Thus, while we do not in general exclude κ d = −1 we can, in the 2HDM-II that we examine here, rule out the arrangement κ u = −κ d = 1 (given by Eq. (5.3)) at 2.7σ from Fig. 2.
More generally, we observe that when κ d, = −1, small deviations away from κ u = 1 have a significant impact, causing cos(β − α) to be large, which is incompatible with the signal strength data. The debate in the literature over the possibility of the wrong-sign limit will hopefully be resolved by future collider data, with suggestions on fruitful search channels made in Refs. [50] and [51], with the latter also including the possibility of a hidden sequential fourth generation of fermions.

Flavour Observables
In this section, we present our fits of the 2HDM-II for selected flavour observables. We show fits for groups of observables individually, and then in Section 7, we combine these with all other observables in the global fit. For each fit, we use the standard (tan β − m H + ) parameter space in a log-log scale. All fits have been performed using the flavio package as described in Section 1.

Leptonic and Semi-leptonic Tree-Level Decays
The Yukawa term in the 2HDM-II Lagrangian (2.17) leads to additional contributions to the tree-level flavour-changing charged transitions. Integrating out the heavy charged Higgs boson H ± , one obtains the effective Hamiltonian describing the d → u −ν or u → d + ν transition in the 2HDM: where the new effective operators are The corresponding Wilson coefficients are related to the 2HDM-II parameters as follows: For use in flavio, we convert the latter coefficients to C du ν SL , C du ν SR , respectively, in the flavio WET basis, finding where the negative sign in front of these coefficients is due to the convention for the effective operators adopted in flavio's basis. Feynman diagrams describing the leptonic and semi-leptonic transitions in the 2HDM are shown in Fig. 4. Figure 4: Diagrams contributing to leptonic (left) and semi-leptonic (right) decays in the 2HDM.
The full list of leptonic and semi-leptonic modes included in our fit is given in Table 6, where the corresponding SM predictions (produced in flavio) are based on [52][53][54][55][56][57][58][59][60]. In addition to previously measured semi-leptonic decay channels, we consider for the first time the B s → D ( * ) s µν µ modes measured recently by the LHCb collaboration [61], and the corresponding B s → D ( * ) s form factors are determined using Lattice QCD [62]. Moreover, we consider the Lepton-Flavour Universality (LFU) The experimental measurements of the latter were found to be in tension with the corresponding SM predictions, giving hints of possible LFU violation. In the 2HDM, this violation is caused by different couplings of the charged Higgs boson to the lepton pair, which are proportional to the lepton mass m . Using the HFLAV averages for R(D ( * ) ) [63], we present in Fig. 5 the allowed (tan β − m H + ) regions at 1 and 2σ levels, where the left (right) plot corresponds to R(D) (R(D * )). Note that the semi-leptonic b → c ν transitions receive negative corrections compared to the SM contributions, therefore R(D) and R(D * ) move further away from the measurements (apart from a narrow region of the (tan β−m H + ) parameter space (see Fig. 5), where the 2HDM-II correction becomes more than the twice the size of the SM contribution proportional to the scalar form factor). Since R(D) is in fact only 1.2σ away from the SM prediction (see Table 6), one finds a wide allowed range of tan β and m H + at the 2σ level (the left plot in Fig. 5). On the other side, R(D * ) is in greater tension (2.8σ) with the SM, therefore one obtains just a narrow region with large tan β and small m H + (the right plot in Fig. 5). The experimental combination of both R(D) and R(D * ) [63] is more than three standard deviations away from the SM prediction, and within 2σ one finds in the 2HDM-II only a very narrow region at very low masses of the charged Higgs, m H + ∼ 1 GeV, which is far from the physical domain. We find that the 2HDM-II is not able to accommodate the experimental data on both R(D) and R(D * ) within 3.5 σ; the corresponding tension in the SM (using flavio) is 3.2 σ.
Combining all the leptonic and semi-leptonic tree-level decay channels indicated in Table 6 yields the allowed range of tan β and m H + at 1 and 2σ levels as shown in Fig. 6.  and π mesons and the hadronic decays of τ leptons to K and π mesons with a tau neutrino as well as R(D) and R(D * ), see Table 6. The lighter contour indicates the allowed parameter space at 2σ confidence level while the darker contour corresponds to 1σ.
In addition, we would like to make an important observation regarding the extraction of the CKM matrix elements from experimental data. Consider e.g. the branching fraction of a leptonic meson decay in the 2HDM, where δ 2HDM is the 2HDM correction factor. Conventionally, the corresponding CKM element is determined from the experimental measurement, assuming B exp = B SM . But if the 2HDM is realistic then B exp = B 2HDM and the CKM element extracted from measurement is actually proportional to (1 + δ 2HDM ). Therefore, the fraction 1/(1 + δ 2HDM ) is then interpreted as the modification factor that the measured CKM element must receive in the 2HDM to be the true CKM element. A similar argument holds for the semi-leptonic meson decays.
Using Eq. (6.5), and a similar relation for the semi-leptonic decays, the size of the modification factor 1/(1 + δ 2HDM ) can be scanned across the 2HDM-II parameter space (tan β − m H + ) to test the significance of the factor for each of the CKM matrix elements, where the accepted values are commonly taken from leptonic and semi-leptonic tree-level decays. We perform tests of this significance for V us , V ub and V cb using the tree-level parameterisation to consider the propagated effects to all CKM elements, and we also take into account unitarity of the CKM matrix. We find negligible effects to all CKM elements for most of the parameter space considered here, and the small areas where more noticeable effects arise are excluded from theoretical constraints, direct searches, and by the global fit in Section 7.
Were the 2HDM-II proven to be physical, it could be the case that there is some small modification to CKM elements that, although is not enough to impact results here, would need to be considered. This would be dependent on the specific values found for m H + and tan β from experiment.

Neutral B-Meson Mixing
Neutral B-meson mixing plays an important role in flavour physics. For a recent discussion on Bmixing in the SM and within New Physics models see, for example, Ref. [64]. In the 2HDM, B-mixing gains additional contributions from new box diagrams where one or both of the virtual W ± bosons are replaced by charged Higgses H ± , see for ∆B = 2 processes can be expressed as [65]: where the operators are defined (for q = d, s) as with α and β denoting the colour indices. The contributions of these operators to B-meson mixing in the 2HDM are well-known throughout literature, see e.g. Refs. [65][66][67]. In our work, for the Wilson coefficients C ( ) k in Eq. (6.6) we use expressions calculated in Ref. [65]. We convert these to the flavio WET basis [12], where we find . In our analysis we will use the averages presented in Ref. [64], which are based on HQET Sum Rule evaluations [68][69][70] and lattice simulations [71][72][73]. The perturbative SM corrections are known and implemented to NLO-QCD accuracy [74]. In the numerical analysis we further use the most recent experimental averages for ∆m d and ∆m s [63] 8 (see Table 6). The results of the fit are shown in Fig. 8, indicating that the mass differences constrain the allowed tan β region from below.

Loop-Level b → s, d Transitions
The effective Hamiltonian for flavour-changing neutral current (FCNC) b → q + + and b → sγ (for q = d, s) processes is defined as [65] H b→q The operators in this Hamiltonian are (6.10) The expressions for 2HDM-II contributions to each operator's Wilson coefficient are given e.g. in Section 3 of [65]. Since all the operators in Eq. (6.9) are defined in the same way as in the flavio WET basis (for either = e, µ), no Wilson coefficient conversions are required here. Some example diagrams describing the b → q + − transitions in the 2HDM are shown in Fig. 9.   leads to a stronger constraint on the charged Higgs mass: m H + ≥ 800 GeV (95% CL) [76]. The experimental average from [63] is based on [79][80][81]. The operators which contribute to the radiative decaysB → X s γ are O 7 , O 8 , given by Eq. (6.10). As expected, our fit for this process provides an important lower bound for the charged Higgs mass, predicted here as (see also Fig. 11) We slightly differ from Ref. [76] in that our 2HDM-II contributions are taken at NLO [82], while they use NNLO results [83] and we also use a different statistical treatment. In addition, as can be seen in Fig. 11, the bound for m H ± becomes even stronger for lower values of tan β.

Leptonic
The FCNC leptonic meson decays B d,s → µ + µ − are particularly sensitive to scalar operator contributions in NP models, and therefore can be excellent probes of the effects of the 2HDM. These processes are particularly worthy of consideration now, due to the recently-produced experimental combinations which we use in our fit [84], and recent ATLAS, CMS, and LHCb results [85][86][87][88][89][90]. The SM prediction is based on a perturbative element [91][92][93] and a non-perturbative determination of the decay constants, see e.g. Ref. [94][95][96]. It has been common in the past to study these decays in the large tan β limit (tan β m t /m b , see e.g. Ref. [97]) where the Yukawa coupling to b quarks is large and there can be simplifications to the 2HDM contributions of the operators which affect this process: O P . However, a big part of the allowed parameter space of the 2HDM-II lies not within the large tan β, thus we use only the general expressions given in Ref. [65]. These generalised expressions bring three further parameters of the 2HDM into play: α, m H 0 , and m A 0 . In Section 7, we will combine all observables in a fit of the five independent parameters: tan β, m H + , m H 0 , m A 0 and cos(β − α). Here, in order to present a 2D fit in the (tan β − m H + ) parameter space, we fix the other three parameters as cos(β − α) = 0, as the alignment limit is preferred (see Section 7) by the global fit, and m H 0 = m A 0 = m H + , which is justified from theoretical constraints (see Section 3). The corresponding contour plot is shown in Fig. 12. From there one can see that the mass of charged Higgs is also constrained from below, roughly: however this bound is not as strong as that from the B → X s γ decay. Furthermore, B d,s → µ + µ − yields a strong correlation between possible values of tan β and m H + .

Semi-leptonic b → s + − Transitions
For the semi-leptonic b → s + − processes we are considering, all of the operators from Eq. (6.10) contribute. Therefore all the 2HDM parameters, m H + , m H 0 , m A 0 , tan β and cos(β − α), will affect these processes.
The lighter contour indicates the allowed parameter space at 2σ confidence level while the darker contour corresponds to 1σ.
The first group of observables are the LFU ratios R K and R K * [98], defined as where the branching fractions are integrated over bins of squared dilepton invariant mass q 2 . These quantities are theoretically very clean, since almost all hadronic effects drop out in the ratio. In the SM, the values of R K ( * ) are very close to 1 with tiny uncertainties. The electromagnetic corrections for R K ( * ) worked out in Ref. [99] were found to be very small, ≈ 1 − 2%. In our analysis, we consider the ten binned LFU ratios listed in Table 8. The recent measurement by LHCb shows a deviation from the SM of 3.1σ for R K + [100]. In Fig. 13 we present the resulting (tan β − m H + ) contour plot for R K ( * ) observables, where we fixed the three other 2HDM parameters at their best fit values for the 10 R K ( * ) bins: cos(β − α) = 0.001, m H 0 = 250 GeV, m A 0 = 9.2 GeV (see Table 2). As one can see, R K ( * ) could only be accommodated within the 2HDM-II only at very low m H + 1 GeV and very high tan β 300, which are beyond the physical domain (see Section 3). Therefore, the 2HDM-II (as well as the SM) is unable to accommodate the current experimental values of R K ( * ) . The contours in Fig. 13 are not only found in regions which are in disagreement with the allowed contours in Figs. 6, 8, 11, 12, but also with the constraints from direct searches [32]. Moreover, in deriving the formulae for the relevant Wilson coefficients [65], it was assumed that m H + is at least of the order of the electroweak scale, meaning that the validity of contours much lower than log 10 [m H + /GeV] ∼ 2 must be questioned. The combined fit of the 10 R K ( * ) observables is in disagreement with our combined fits of all observables to 4.2 σ confidence. Due to this disagreement between R K ( * ) and other flavour observables,   Table 2. The lighter contour indicates the allowed parameter space at 2σ confidence level while the darker contour corresponds to 1σ.
we work with two approaches for our fits: excluding and including R K ( * ) .

Results
Combining all the observables collected in Table 4 (Higgs signal strengths), Table 5 (S, T, U ), Table 6 (flavour observables without b → s + − ), Table 7 (b → s + − observables) and Table 8 (R K ( * ) ) one obtains the results shown in Table 2. We emphasize that in searching for the best-fit point of the corresponding likelihood function we have applied the theory constraints for the Higgs mass differences (see Table 1). As one can see from Table 2, including all the observables yields a relatively poor fit with a low p ≈ 1.5%-value. This is due to the fact that the R K ( * ) and R(D ( * ) ) quantities cannot be accommodated within the 2HDM-II, as was already discussed in Sections 6.1 and 6.3.3. Nevertheless the 2HDM-II still achieves a higher level of agreement than the SM and we find a corresponding pull of 2.3σ. Not considering the observables R K ( * ) yields a better fit for the 2HDM-II  Table 7 excluding the R K ( * ) observables and fixing the additional parameters to the best fit point from the b → s + − fit in Table 2. The lighter contour indicates the allowed parameter space at 2σ confidence level while the darker contour corresponds to 1σ.
with a significantly higher p ≈ 7%-value, again marking a better performance than the SM, with a corresponding pull of 1.8σ.
This bound is considerably stronger than the one we obtained in Eq. (5.4) using only data on Higgs signal strengths. Regarding the wrong sign limit solution, we found using only the Higgs signal strengths this to be allowed at 2.7σ or above; using the global fit of all observables (excluding R K ( * ) ), this is now excluded within 6σ of our best fit point with p-value ≈ 7%. We also perform fits to all observables requiring the wrong sign limit as defined in Eq. (5.3), where one can see the quality of these fits is significantly worsened compared to the general global fit with no input constraint on cos(β − α), yielding a p-value ≈ 0.1%.

The best fit further favours values of m H
TeV and tan β ≈ 4. Performing a scan around the best-fit point (excluding R K ( * ) ) and keeping in mind the theory constraints in From the 2σ-level onwards we do not see an upper bound within the area considered (m H + ≤ 100 TeV). The parameter tan β has a strong correlation with m H + . For the best fit values of cos(β − α), m H 0 and m A 0 , we find the allowed regions in the (tan β − m H + ) plane given in Fig. 16 and Fig. 17 (the latter zoomed into the best-fit point and shown in combination with the theory constraints). Assuming both the alignment limit, cos(β − α) = 0, and the degeneracy for the Higgs masses, m H + = m H 0 = m A 0 , the contour plot of Fig. 18 shows that tan β is strongly constrained for low m H + ∼ 1 TeV while the bounds on tan β become weaker with an increase of m H + up to 100 TeV. We stress again that these constraints above are obtained from the fit by excluding R K ( * ) . The inclusion of R K ( * ) would not significantly change the allowed parameter regions, while it would considerably worsen the quality of the fit. Furthermore, we investigate several constrained scenarios in the fit, as indicated in Table 2. First, we checked that the wrong-sign limit is absolutely disfavored by the global fit. Second, in the exact alignment limit the fit yields similar results as in the general case. Finally, assuming both the alignment limit, cos(β − α) = 0, and the degeneracy for the Higgs masses, m H + = m H 0 = m A 0 , again leads to similar results due to the preference of both of these assumptions by the global fit and theory constraints.

Comment on the Anomalous Magnetic Moment of the Muon
The anomalous magnetic moment of the muon, measures the deviation of the muon gyromagnetic moment g µ from 2 due to quantum loop effects. This can be accurately measured in experiment, with recent measurements at Fermilab [13] confirming the older BNL value [14]. It can also be accurately predicted in the Standard Model and the Theory Initiative (WP [15]) yields a 4.2σ discrepancy with the updated combined experimental value. We will also consider a recent Lattice QCD evaluation (BMW Collaboration [16]) that finds the g µ − 2 discrepancy to be only 1.6σ. This deviation in a µ can be a strong motivation for potential new physics signals, such as supersymmetric particle loops or other BSM contributions, see e.g. Ref. [116]. 2HDM contributions to a µ form a subset of SUSY contributions to the anomalous magnetic moment and they have been recently studied e.g. in Refs. [117][118][119].
Matching flavio's WET-3 basis [12], we write the effective Hamiltonian for → γ processes as  Table 2: Best fit points of 2HDM-II parameter fits for various groups of observables using the constraints from theory to inform the physical parameter values.
The 2HDM contribution to a µ is then defined as where we take the contributions at one-and two-loop level given in [120]. The two-loop contributions come from the Barr-Zee diagrams ( [121], see Fig. 19) with a fermion or boson loop which can give rise to significant contributions to a 2HDM µ . All Higgs bosons are present in these loops, therefore a 2HDM µ depends not only on tan β and m H + , but also on m A 0 and m H 0 , with a further dependence on the mixing angle α. From Fig. 20 we find that the 2HDM-II cannot explain the discrepancy between the experimental number for the anomalous magnetic moment of the muon and the white paper prediction as the required values of tan β lie in the non-perturbative region. Obviously the BMW prediction yields a consistency of experiment and theory within 2σ and thus a large portion of the 2HDM-II parameter space is still allowed.  Table 2. The contours indicate the allowed parameter space at 1, 2, 3, 4, 5σ confidence level going from darker to lighter.   Table 2. The contours indicate the allowed parameter space at 1, 2, 3, 4, 5σ confidence level going from darker to lighter.

Electroweak Phase Transition
The electroweak phase transition (EWPT) is a mechanism which can generate the baryon asymmetry we observe in the Universe today through the process of EW baryogenesis, with the requirement that  Table 2 and including the restrictions from the theory constraints in Table 1 (shown by the red dashed lines). The contours indicate the allowed parameter space at 1, 2, 3, 4, 5σ confidence level going from darker to lighter.
The contours indicate the allowed parameter space at 1, 2, 3, 4, 5σ confidence level going from darker to lighter.  the EWPT is of strong first order (SFO). In the SM this could be achieved for m h 0 70 GeV [122], however the measurement of the SM Higgs boson mass as m h 0 = 125.1 ± 0.14 GeV [21] means we now require BSM physics to achieve a SFOEWPT; a 2HDM is in principle capable of generating this. For recent work testing the 2HDM for a SFOEWPT across large regions in its parameter space, see e.g. Refs. [19,123,124], in which the authors find regions where the 2HDM-II could support a SFOEWPT, however only for Higgs masses less than 1 TeV. In [125], the strength of the EWPT in the 2HDM-II is tested for Higgs masses 600 GeV; a FOEWPT can be found for these masses but it is not strong enough to support EW baryogenesis. Furthermore, a SFOEWPT appears to favour a mass splitting (∼ 200 GeV [123]) between H + and H 0 , and m A 0 ∼ m H + [123,124]. This could still be allowed by our results in Section 3, however only for low charged Higgs masses, in tension with the global fit of Fig. 16.
In this work, we use the BSMPT package [17,18] to calculate the strength of the EWPT for selected parameter points across our fit. The BSMPT package calculates the strength of the EWPT for the 2HDM in the lambda parameter basis (λ i , m 2 12 , tan β). We have so far primarily worked in the mass basis (m A 0 , m H 0 , m h 0 , m H + , tan β, cos(β − α)), and so must convert to the lambda basis. In the alignment limit, a sufficient but not necessary condition for vacuum stability is [22] We then use this, along with Eqs. (2.11)-(2.15), to perform the conversion, while making the assumption that cos(β − α) is small, which is of course in line with our constraints. This is only appropriate when the conditions of Eq. (7.6) are fulfilled, i.e. primarily for m H 0 m A 0 , m H + . The BSMPT package parameterises the strength of the EWPT in the 2HDM by calculating where ω c is the high-temperature VEV (ω 2 c = ω 2 1 (T c ) + ω 2 2 (T c ) and ω i (T = 0) = v i ) at the critical temperature T c . The BSMPT package only passes a result if ξ c > 0, with the desired SFOEWPT indicated by ξ c > 1.
We choose M = 50 TeV as a benchmark for more extreme masses, the best-fit point of our degenerate mass fit, and lower limits in m H + at ∼ 2, 3, 4, 5σ where we select more favourable mass splittings as points to test for a SFOEWPT. We also show some variations in tan β at the best-fitting degenerate masses. In Table 3, we present our selected parameter points in both bases, and then also the results produced by the BSMPT package. As is generally favoured by our fits, for higher masses we mostly test parameter points in the limit of degenerate masses, since this will satisfy the conditions of Eq. (7.6) and the small mass splittings allowed at higher mass scales have only small effects, as is shown in one benchmark point. However for lower masses where these splittings are more significant, we consider further benchmark points allowing for mass splitting more favourable for a SFOEWPT [123], and find that demanding degeneracy for these lower masses is a poorer choice.
We find that for most of our allowed parameter space, a SFOEWPT is not possible due to the large mass scales involved. From the benchmark points at higher masses, we find a maximum of ξ c = 0.17; much smaller than required. For our benchmark points at lower masses however, the outlook improves, where we first find a ξ c ≥ 1 for a parameter point at ∼ 2σ from the global best fit point. On the other hand, we note that all points in Table 3 with a ratio ξ c ≥ 1 require some of the couplings |λ i | > 4, which could be in conflict with perturbativity, depending on what constraint we apply; see the beginning of Section 3.1. As we have only chosen a few benchmark points, we cannot rule out the possibility of a ξ c ≥ 1 at some point closer to our best fit, however this seems unlikely given the trends shown here and in other works. Therefore to accommodate a SFOEWPT in the 2HDM-II, we would require much lower Higgs masses than are favoured by the global fit. But for such regions of the parameter space the 2HDM-II will no longer perform better in the overall fit compared to the SM.   Table 3: Table of results for the EWPT in the 2HDM. We consider the exact alignment limit, cos(β − α) = 0, and in addition we employ the condition by Eq. (7.6). In this limit, as follows from Eqs. (2.11)-(2.15), one gets: λ 1 = λ 2 = m 2 h 0 /v 2 = 0.26, i.e. λ 1,2 are fixed and independent of the Higgs masses and tan β (and therefore not again shown in the table). In addition, one obtains that λ 3,4,5 are independent of tan β, and λ 3 + λ 4 + λ 5 = m 2 h 0 /v 2 = 0.26. For similar discussions of the lambda basis, see e.g. [126].

Conclusion
In this work we have investigated theoretical bounds on the 2HDM-II stemming from perturbativity, unitarity and vacuum stability as well as experimental constraints stemming from electroweak precision, Higgs signal strengths, flavour observables and the anomalous magnetic moment of the muon. All in all, we find that the 2HDM-II can accommodate the data better than the SM. The best fit point for the 2HDM-II lies around For the charged Higgs mass we find a lower limit of 680 GeV at 3σ and the remaining Higgses have to be largely degenerate -this requirement becomes stronger the heavier the Higgs masses become. For low values of the charged Higgs mass, tan β is severely constrained to lie around a value of 4, while for high values of the charged Higgs mass (around 50 TeV) tan β should be in the region between 0.1 and 300. We further find that the wrong-sign limit is disfavoured by Higgs signal strengths and excluded by the global fit by more than 5σ and that only small deviations from the alignment limit are possible.
Experiments currently observe some deviations from SM predictions, in particular in R(D ( * ) ), R K ( * ) and the anomalous magnetic moment of the muon. The tensions in R(D ( * ) ) will either be enhanced in the 2HDM-II or to be resolved, they would require parameter regions in the (tan β−m H + ) plane which are excluded by other flavour observables. Eliminating the tension in R K ( * ) would require charged Higgs masses below 1 GeV and non-perturbative values of tan β. The 2HDM-II can also not resolve the discrepancy between the measurement of the magnetic moment of the muon and the white paper prediction.
Finally we found that the allowed parameter space of the 2HDM-II strongly disfavours the existence of a sufficiently strong first order phase transition for EW baryogenesis. This requires mass splittings discouraged by theoretical constraints, and low masses of the charged Higgs -at least 2σ away from our best fit point, which makes this model theoretically much less attractive.

A Inputs and Observables
We start from Eq. (2.10) and rearrange to give an expression for m 2 12 . Using small angle approximations on δ = β − α − π/2, we expand to second order, and find that where a = 2s 2β ∆t β Γ −2 β + 2c 2β ∆t β Γ −1 β − 2s 2β ∆t β Γ 0 β − s 2βλ1 c 2 β + s 2βλ2 s 2 having defined ∆t β = tan(β) − 1 tan(β) . (B.10) By trigonometric manipulation we find that c = 0. We are therefore left with the two solutions where the first solution corresponds to the alignment limit. Note that the second solution is only valid for small δ as a second-order Taylor expansion was used in its derivation.