The Flavourful Present and Future of 2HDMs at the Collider Energy Frontier

We study the intersection of flavour and collider physics for Two-Higgs-Doublet models of Type I and II. Drawing from the flavour precision-LHC exotics search complementarity, we also provide a projection of the future sensitivity that can be achieved in light of currently available analyses. On the one hand, we find that the parameter space of the 2HDM can be explored significantly further with more data from the LHC with some complementarity with flavour physics. On the other hand, flavour physics results alongside their projections remain powerful tools to constrain the model space in regions where direct sensitivity to new states via exotics searches is lost. Our results further highlight the recently observed flavour physics anomalies as important drivers of new physics searches in the future; we also touch on implications for a strong first order electroweak phase transition.


Introduction
Searches for physics beyond the Standard Model (BSM) are key pillars of the current particle physics programme that stretches across many different experimental arenas. Shortcomings of the Standard Model (SM), e.g. the possible lack of enough CP violation or the insufficiently strong electroweak phase transition to address the Sakharov criteria [1], provide the motivation for extending the interaction and particle spectrum of the SM. Especially well-motivated scenarios along these lines are Higgs sector extensions. While the existence of a non-trivial electroweak vacuum structure had been verified with the discovery of the W and Z bosons at UA1/UA2, it was only in 2012 that the discovery of the Higgs boson [2,3] confirmed the mechanism of spontaneous symmetry breaking [4][5][6] through the Higgs boson's couplings to the massive gauge bosons. Subsequent investigations of the 125 GeV Higgs boson after its discovery have shown that this new state closely follows the SM expectation. Nonetheless, current experimental constraints only limit modification of the Higgs boson's interactions to be below ∼ 10% [7,8] so there are reasonable margins for new, weakly-coupled modifications of the Higgs sector. Consistency with electroweak fits [9] provides an additional constraint for extensions of the symmetry breaking Higgs sector that can be competitive with current direct measurements.
Out of the Higgs sector extensions that are typically considered for BSM investigations, two Higgs doublet models (2HDMs) are particularly well-motivated theories (for reviews see Refs. [10,11]). On the one hand, they implement custodial isospin analogous to the SM Higgs field, thus avoiding tensions and fine-tuning with electroweak precision measurements that occur in higher dimensional representations of electroweak SU (2) L ×U (1) Y breaking [12]. On the other hand, 2HDMs can be considered as harbingers of supersymmetric theories, which demand a second Higgs doublet due to holomorphy, and (non-)perturbative anomaly cancellation to maximise the S-matrix symmetry.
Leaving aside the theoretical motivation for introducing a second weak doublet contributing to electroweak symmetry breaking and fermion mass generation, the experimental avenues to constrain or even observe a 2HDM extension of the SM are plentiful: links of fermion vs. gauge boson mass generation provide a motivation to combine precision flavour physics investigations with those performed at the high energy frontier, chiefly at the Large Hadron Collider (LHC). In fact, a recent investigation into the flavour and Higgs signal strength constraints of the 2HDM [13] demonstrated the power of flavour measurements in driving BSM mass scales into a region where the direct detection at the LHC could become challenging. Additionally, current flavour anomalies highlight flavour physics as a particularly relevant area to inform future investigations at the high energy frontier.
Building on Ref. [13], we investigate the intersection of flavour and collider physics for the 2HDMs of type I and II in this work, with a view towards the High Luminosity (HL) LHC phase. Clarifying and extrapolating the currently available search strategies in flavour physics and for the exotic states that are predicted in the 2HDM enables us to discuss flavour-collider complementarity and identify regions that require targeting to enhance the BSM potential during the upcoming LHC runs.
We organise this work as follows: in Section 2, we review the basics of the 2HDM to make this work self-consistent. Section 3 provides details on theoretical and electroweak constraints that set the baseline of our investigation, then in Section 4 we discuss the status of the SM Higgs signal strengths in the 2HDM type I (2HDM-I). In Section 5 we discuss in detail the flavour constraints on the 2HDM-I, extending the work of Ref. [13], before we move on to contrast LHC searches and measurements in Section 6 for the 2HDM-I and 2HDM-II. We also comment on expected improvements of both flavour and collider constraints before discussing cosmological implications in Section 7. We provide a summary and conclusions in Section 8.

The Two Higgs Doublet Model
We follow the notation of Refs. [10,11], where, instead of using the SM Higgs doublet and its (symplectic) charge-conjugated version to gain mass terms for both up-type and down-type quarks, one can introduce two distinct doublets with i = 1, 2. After electroweak symmetry breaking (EWSB), these doublets lead to 5 physical Higgs bosons (the other 3 degrees of freedom construct the longitudinal W ± , Z 0 boson polarisations). One finds two charged Higgs, H ± , two neutral scalars, h 0 , H 0 , and a neutral pseudoscalar, A 0 . We take h 0 to be the lighter of the two neutral scalars and the currently observed Higgs particle with mass 125.25 ± 0.17 GeV [14]. The potential for a general 2HDM with a softly broken Z 2 symmetry is [10] We can use this potential in its "lambda basis" to construct mass terms for each of the physical Higgs particles, which, alongside the vacuum expectation value (VEV) v 2 = v 2 1 + v 2 2 , the ratio of VEVs tan β = v 2 /v 1 , the mixing parameter cos(β −α), and the softly Z 2 breaking term m 12 form the "mass basis" which more simply translates to physical observables, see e.g. Refs. [13,[15][16][17][18] for transformations between bases.
The Yukawa couplings in the 2HDM can be generally expressed as In general the 2HDM Yukawa couplings need not be flavour diagonal and so it is possible to allow tree-level flavour-changing neutral currents (FCNCs) which are absent in the SM. In order to maintain this natural flavour conservation, one can consider fermions coupling to the Higgs doublets in specific ways. The four types of 2HDM with natural flavour conservation are described in Table 1. After spontaneous symmetry breaking, we can write the Yukawa sector of the 2HDM Lagrangian as (see e.g. Refs. [10,19]) Imposing the fermion interactions to specific doublets in Table 1, the ξ coupling strengths are expressed in Table 2 Table 2. Coupling strengths ξ in each type of 2HDM between the Higgs particles and fermions.

Perturbativity, Unitarity, and Vacuum Stability
Perturbativity in the scalar sector can be simply expressed as [20,21] As in Ref. [13] we also consider the less conservative bounds of |λ i | ≤ 4 -inspired by the results in Refs. [22,23]. Looking at the couplings in Eq. (2.4), one can further consider the perturbativity constraints from the Yukawa sector: where we 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. For the upper bound we adopt the same tan β = 10 +2.5 ≈ 300 as in Ref. [13], for ease of comparison of results between 2HDM-II [13] and 2HDM-I (this work). We apply the conditions for a stable vacuum as set out in Ref. [24], whilst demanding the vacuum to be the global minimum of the potential [25]. We also consider the conditions from tree-level unitarity, see Refs. [21,26] 2 , alongside NLO unitarity and the condition that NLO corrections to partial wave amplitudes are suppressed relative to LO contributions, see Refs. [22,23].
We get the same constraints on the 2HDM parameters as found in our previous study from theory, where these are not dependent on the specific couplings defining the different 1 Note that we differ from Ref. [10] by a sign in the H + Yukawa interaction which in turn leads to a difference in the fermion coupling strengths. This is the same convention as found in e.g. Ref. [19]. 2 For a similar discussion in an alternative lambda basis see Ref. [27].
types of 2HDM. We refer to Fig. 1 of Ref. [13] for these results, where the main implication is a close mass degeneracy for the new Higgses enforced from ∼ 1 TeV and becomes stricter as the mass scale increases.

Oblique corrections
The expressions for the electroweak precision observables (EWPOs): S, T , and U [28,29] (also referred to as oblique corrections) in the 2HDM (derived from Ref. [30]) are explicitly given in Appendix C of Ref. [13]. New physics generally has only a small effect on U in comparison to the S and T parameters as the latter correspond to distinct dimension 6 operators in the effective field theory expansion, in contrast to U [31]. It is therefore justified to follow the approach that neglects U from fits by setting U = 0 as outlined in Ref. [9,14]. 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
Many measurements of the properties of the observed Higgs boson have been made. Of these, the most relevant for our discussion are the Higgs signal strengths µ f i , which are defined as a ratio between experimental and SM values of the product of the cross section and branching fraction in a given channel with production mode i and decay products f : Measurements of µ f i essentially serve as an indicator of how closely the observed Higgs matches the SM expectations. We include results from 31 channels [8,[32][33][34][35][36][37] (collected in Table 4 of Ref. [13]).
If the observed Higgs is not that of the SM, but part of an extended Higgs sector such as the 2HDM, its couplings to other particles will differ from the SM, a difference that will feed through into the values of the signal strengths. This is the case in the 2HDM-I, with the couplings being modified by the factors shown in Table 2. It is therefore possible to use the signal strength measurements to constrain the parameters on which these couplings depend, namely the mixing angles. We do this by performing a fit using analytical calculations of the signal strengths in the 2HDM-I as functions of the input parameters tan β and cos(β − α) from expressions given in Ref. [11], with results presented in Fig. 1. Opposed to the 2HDM-II case (see Ref. [13]) we find that more sizeable deviations from the alignment limit cos(β − α) = 0 are allowed for the 2HDM-I. We observe that larger deviations from the alignment limit are possible with increasing tan β, and we find | cos(β − α)| ≤ 0.21 (0.4) at 2 (5)σ for sufficiently large tan β ≥ 10.

Flavour Constraints in 2HDM-I
In Section 5.1, we cover the constraints on the 2HDM-I parameter space from flavour observables, and then in Section 5.2 we combine these with the Higgs signal strengths and the EWPOs. We follow the procedure of the earlier study [13] in the 2HDM-II, where we now exchange the Type II couplings for the Type I. All observables we consider here are listed with SM predictions and experimental measurements in Tables 4-8 of Ref. [13]; the measurements are collected from various results throughout literature [8,14,32,33, and the SM predictions are calculated in flavio unless otherwise stated.

Individual Fits
For the tree-level (semi-)leptonic flavour-changing charged transitions, the 2HDM-I contributes through the effective operators with P L(R) = (1 ∓ γ 5 )/2. The Wilson coefficients of these operators in terms of the 2HDM-I parameters are For a description of all operators gaining contributions from the 2HDM for the observables we consider and the translation of basis for these, we refer to [13], taking expressions for the Wilson coefficients from Refs. [19,67,68] where these are expressed in terms of general 2HDMs and one can simply insert the appropriate couplings for the 2HDM-I. First we consider the Lepton-Flavour Universality (LFU) observables R(D ( * ) ) ≡ B(B → D ( * ) τν τ )/B(B → D ( * ) ν ), where = e, µ. The latter implies a 2.8σ tension between experimental measurements and SM predictions, and the combined tension of the two between experiment and the SM is 3.2σ (using flavio). We show in Fig. 2 the allowed 2HDM-I parameter space for these two observables in individual fits. Within the theoretically-motivated parameter limits of our fits, we find that R(D) allows a large part of our parameter space within 2σ, whereas R(D * ) allows most of our parameter space at 2.8σ or above. In fact, for both observables, we find the best fit for m H + ∼ 1 GeV which is outside the physical domain of our model. Within the considered limits of our parameters, we find a minimum of 3.5σ combined tension between the 2HDM-I and experiment for R(D ( * ) ). In Fig. 3, we present the combined fit of tree-level flavour-changing charged currents including leptonic and semi-leptonic decays of mesons, hadronic decays of τ leptons, and R(D ( * ) ); see Table 6 of Ref. [13] for all channels considered. The SM predictions in flavio for the leptonic and semi-leptonic channels are based on Refs. [69][70][71][72][73][74][75][76][77], and the new 2HDM-I contributions are described by Eq. (5.2).  Next we consider the mixing of neutral B d,s mesons. The mass differences of this meson mixing are experimentally known very precisely at the level ∼ O(0.1%) [14]. On the theory side, however, the uncertainties are still dominated by the non-perturbative determinations of the matrix elements of the ∆B = 2 operators. We use the averages presented in Ref. [78] combining HQET Sum Rules [79][80][81] and lattice calculations [82][83][84], yielding a theory precision of O(5%). Terms proportional to the tH − coupling dominate the 2HDM contributions in the box diagrams with the down-type quark effects largely mass-suppressed. The predictions in the 2HDM-I for B-meson mixing are therefore very similar to those in the 2HDM-II where the dominating up-type H ± couplings remain the same. The combined fit for ∆m d,s is shown in Fig. 4, where we find a correlation between tan β and m H + constraining tan β from below. We take the expressions for the ∆B = 2 Wilson coefficients from Ref. [19], converting to the flavio basis as described in Ref. [13].
A signature observable for new physics analysis is the branching ratio of theB → X s γ decay. In the SM, this is known to NNLO in QCD [85] (based on the previous Refs. [86,87]). We take the experimental average [43] formed from Refs. [88][89][90]. In the 2HDM-II,B → X s γ provides a distinct constraint on the lower bound for the charged Higgs mass m H + ; in the Type I however, this is much more correlated with tan β, and we therefore cannot find a clear constraint from this observable alone. In Fig. 5 we show the fit to B(B → X s γ)| Eγ >1.6 GeV in the 2HDM-I where we calculate the contributions to the Wilson coefficents C 7 , C 8 at NLO using Ref. [68].
We also consider the FCNC processes B d,s → µ + µ − , which are sensitive to BSM contributions to scalar operators, making them important observables to test the 2HDM. In recent years, ATLAS, CMS, and LHCb [91][92][93][94][95][96] have all produced measurements for B s → µ + µ − , and developed further the upper limit on B d → µ + µ − . In our analysis we make use of the combination of these measurements from Ref. [47]. The theory prediction is formed of the perturbative calculations [97][98][99] and the determinations of the non-perturbative  Figure 6. Contour plot of the allowed 2HDM-I parameter space in the (tan β − m H + ) plane for the combined fit to (B d,s → µ + µ − ), taken in the limits of alignment (cos(β − α) = 0) and degenerate masses (m H 0 = m A 0 = m H + ). The darker contour indicates allowed parameter space at 1σ confidence, and the lighter at 2σ. decay constants, for example [100][101][102]. The 2HDM Wilson coefficients (taken again from Ref. [19]) for the operators contributing to these processes, O In this approach, we find here similarly toB → X s γ a strong correlation between tan β and m H + . Later in our global fit (Section 5.2), we will discuss further the dependence on these additional parameters.
The LFU ratios R K and R K * [103] have become well-studied processes in recent years (see e.g. Ref. [104][105][106][107][108][109][110][111]), driven by LHCb measurements [64,65] finding discrepancies with  the SM, most noticeably the deviation of 3.1σ in R K + [63]. Model-independent analyses (for example, Refs. [112][113][114][115]) show these processes favour vector-like new physics contributing to O 9,10 for their deviations to be resolved. In Table 3 we show the best fit point only including the 10 R K ( * ) bins we consider, where we find that this fit favours the limits This result is purely numerics-driven; we note that it lies outside the bounds of our model given by theoretical considerations. In the large tan β limit we recover the SM predictions, as the 2HDM-I induced coupling deviations will go to zero. With the SM result equating to our best fit point, no point in our parameter space can reduce the tensions in these observables beyond their current status. In Fig. 7, we show our combined fit for the 10 R K ( * ) bins within a more similar parameter region to other fits, where we see that most of our parameter space considered lies within 1σ of the best fit (or SM) result. A discovered 2HDM-I within the physical bounds of the model could only increase the tensions in R K ( * ) as we find currently them in the SM.
In Fig. 8, we consider more b → s + − observables (see Table 7 in Ref. [13]), where some of these also find tension with experiment in the SM. For detailed descriptions and analyses focused specifically on this interesting group of observables, see e.g. Refs. [47,72,105,107,[110][111][112][113][116][117][118][119][120][121][122][123][124]. Similarly to R K and R K * above, these processes tend to favour new physics (NP) contributions from O ( ) 9,10 . In the 2HDM-I, we find that the contributions to these operators (and also to O ( ) S,P ) can be sufficient to reduce the tension with experiment for some of these observables. This effect is most significant for m H 0 ∼ m A 0 ∼ m H + ∼ 1000 GeV as suggested by the best-fit point in Table 3, where we find an improvement over the SM with a pull of 2.6σ (excluding R K ( * ) ); we then find a large portion of our parameter space within 1σ confidence of this result. Varying the fit parameters around the best fit point, we also find that this group of observables imposes a lower bound on tan β as tan β 0.14 (0.16) at 2σ (1σ), (5.3) which is compatible with the lower bound coming from perturbativity in Section 3.

Global Fits
We now consider our global fits to the 2HDM-I combining flavour observables, Higgs signal strengths, and the EWPOs S, T, U . Within the contours of a global fit, the LFU ratios R D ( * ) , R K ( * ) still have a tension ∼ 3.5σ with experiment that is worsened compared to the SM. 3 This motivates two scenarios for the global fits: either including these LFU ratios, or considering a fit excluding them. In Table 3, we summarise the results of the various fits performed, indicating their best-fit parameter points and the corresponding statistics. In Fig. 9, we show the fit to all observables (excluding R D ( * ) ,K ( * ) ) where we fix the additional 2HDM parameters to their best fit point as shown in Table 3. In Fig. 10, we show the fit to all observables (excluding R D ( * ) ,K ( * ) ) in the alignment limit and with degenerate masses. We find the shapes of the contours to be much the same for excluding/including R D ( * ) ,K ( * ) , so we do not show both sets of plots here; the impact of including R D ( * ) ,K ( * ) is clearer in Table 3, where we find distinctly poorer quality of fit (low p-value) in this scenario. Excluding R D ( * ) ,K ( * ) , we find a pull of 1.9σ improvement from the SM to 2HDM-I; including R D ( * ) ,K ( * ) , the improvement is 1.5σ. From Table 3, we see that cos(β − α) = 0 (the alignment limit) is closely favoured from the fits, although Fig. 9 shows that larger deviations from alignment are allowed than in the 2HDM-II. We find preference in the global fits for tan β ∼ O(100) and m H + ∼ where a closer mass degeneracy at this scale is enforced from theory constraints. Within the parameter region found we scan around the best fit point to find the full constraints on our parameters. As one can see throughout this section, much of the small tan β parameter space is excluded by most observables individually and globally. We cannot state an explicit bound on tan β from the global fit however as the lower limit is clearly correlated with the charged Higgs mass. For the mass of the charged Higgs, within the parameter region we test, we find bounds only at 1σ confidence: In our global fit of all observables, our constraints on cos(β − α) are much improved from considering the Higgs signal strengths alone, where we find the constraints in the full analysis to follow very closely to the representative plot in Fig. 9: | cos(β − α)| ≤ 0.21 at 5σ.

Comment on the Muon Anomalous Magnetic Moment
It is worthwhile to comment on the anomalous magnetic moment of the muon, a µ , which we do not include in our overall fit. The recent results from Run 1 at Fermilab [125] confirmed the previous result from BNL [126], where the combined experimental value now yields a 4.2σ deviation from the SM result predicted by the Theory Initiative White Paper (WP) [127] (based on Refs. [128][129][130][131][132][133][134][135][136][137][138][139][140][141][142][143][144][145][146][147]). There is also a competing Lattice QCD prediction for the a HVP µ contribution from the BMW Collaboration [148] which results in only a 1.6σ discrepancy from experiment. 4 Here, we fit a µ to the 2HDM-I parameters considering both the WP and the BMW SM scenarios, shown in Fig. 11. We take the one-and two-loop contributions of the 2HDM to a µ from Ref. [67] where one can insert the expressions for the 2HDM-I couplings and we then convert the results to flavio's WET-3 basis as described in Section 7.2 of Ref. [13]. As before, to present information in terms of two-dimensional fits, we fix the additional 2HDM parameters cos(β − α) = 0 and m H + = m H 0 = m A 0 as motivated by theory. We find that in the 2HDM-I, a µ strongly favours low m H + ∼ 100 MeV, low tan β ∼ 0.5, where the tension between experiment and theory using either SM prediction can be reduced to much less than 1σ. This result is however below the physical domain for m H + . Within the physical domain, we find that m H + 100 GeV can yield a tension between the WP and experiment of less than 5σ, with a minimum of 4.3σ, and between BMW and experiment less than a 2σ tension, with a minimum of 1.6σ. In the left plot, the SM prediction taken from the theory initiative is used; in the right, the SM result from the BMW collaboration. Contours are plotted representing allowed parameter space at 1,2,3,4,5σ confidence from darkest to lightest: in the left plot, only the 5σ contour is visible; in the right, the 2, 3, 4, 5σ contours.

2HDM-I and II Flavour Prospects at Future Colliders
Finally, we would like to discuss the prospects of the 2HDM flavour sector at future colliders. An approximate prediction for B-physics at future colliders such as the HL-LHC or beyond is that the precision of the measurements will roughly double from that of today's [158,159] We perform this process on ∼200 of our 275 observables which are identified as being "experimentally-limited" in their χ 2 contribution: that is, their present χ 2 contribution is dominated by the size of their experimental uncertainty. It is assumed in these tests that the theoretical calculations and precision do not change significantly by the time of these future measurements.
In Table 4, we list the comparison of the χ 2 best fits in the SM, 2HDM-I, 2HDM-II for first the present situation and then each future scenario considered. We also list the lower bound of m H + fromB → X s γ in the 2HDM-II for each scenario, as this has been an important bound in the history of the 2HDM and can remain so in the future. It is unlikely that any one of the exact scenarios would come to pass from future measurements for such a large group of observables, however comparison between all 6 scenarios sheds some light on the future suitability of the SM and the 2HDMs. In the present, we find that both the 2HDM-I and 2HDM-II perform better than the SM by ∼ 2σ. In all future scenarios, the 2HDM-II still performs better than the SM, with a pull > 1σ in all but one of these. The 2HDM-I's performance over the SM is narrowed significantly however, where the best future scenarios for this have a reduced pull of 1σ. In the "Ideal" future, the SM, 2HDM-I, and 2HDM-II all significantly improve upon their current positions; the majority of observables in this "Ideal" scenario simply move much closer to their SM predictions, where the small contributions from the 2HDMs either are not significant or are sufficient to improve upon the remaining small tensions with the SM.
Overall, it is found that increased precision of flavour observables at future colliders can lead to much poorer fits in both the SM and the 2HDM. The exception to this is when measurements are shifted much closer to the SM predictions and the small 2HDM contributions can be sufficient to resolve much of the remaining differences. While the future measurements' central values may not follow the scenarios considered here, the limiting factor in the χ 2 values in Table 4 is the increased precision of measurements. Unless future measurements would indeed conform closely to our "Ideal" future scenario, the performances of the SM, 2HDM-I, 2HDM-II would all become very poor. While assuming no significant gaps in our theory predictions, this suggests some New Physics other than the 2HDM should be explored.  Table 4. Summary of future prediction scenarios and their χ 2 best fits. For the 2HDMs in brackets are the pulls from the SM where a positive value indicates an improvement and a negative value a worsening. The results are shown for excluding R D ( * ) ,K ( * ) . Also shown is the variation of the lower bound of the charged Higgs mass in the 2HDM-II fromB → X s γ.

LHC Constraints, Flavour Comparison and Outlook
the years since the Higgs discovery, many precise measurements of the 125 GeV scalar have been made (see above), while searches for other, more exotic, Higgs bosons have continued alongside the precision measurements of the observed Higgs (see also Refs. [160,161]). These precision measurements closely match the phenomenology predicted of the SM Higgs. In this Section we look to combine the wealth of experimental data on direct searches for new Higgs bosons to further constrain the allowed parameter space of the 2HDM, independently from the constraints in the prior sections, before combining these results in Section 6.3. We also extrapolate the current LHC bounds to the expected results at the HL-LHC operating at 3 ab −1 in the same channels. As in the rest of this work, the parameters of primary interest here are the masses of the new Higgs bosons and the angles tan β and cos(β − α). Following the results of Section 5.2, we examine the favoured scenario in which the new Higgs bosons have degenerate masses and the alignment limit holds.

Current Collider Bounds
In order to leverage the power of the vast array of search data we use the packages HiggsBounds [162][163][164][165][166][167][168] and 2HDecay [169][170][171][172][173][174]. We use 2HDecay to calculate the branching ratios and decay widths of each 2HDM Higgs boson for a given point in parameter space, and interface this data with HiggsBounds to check if the search data excludes such a point. For the key channel of H + production in association with tb, we use MadGraph5 aMC@NLO [175] to generate cross sections which also form part of the input to HiggsBounds, along with the couplings of the Higgs bosons as given in Table 2. We perform scans using 50000 randomly generated points in the parameter space following this method and show the allowed points for both type I and II 2HDMs in Fig. 12. The historic baseline sensitivity from LEP [176] excludes charged Higgs masses below 72.5 GeV and 80 GeV for the Type I and Type II 2HDM models respectively. In the 2HDM-I scan on the left of Fig. 12, at low values of tan β, leptonic decays of the neutral Higgses (e.g. H 0 → µ + µ − in Ref. [177]) exclude lower masses. Sensitivity to H + → τ + ν τ peaks at ∼ 85 GeV [178] and then the H + → tb [179] decay provides the exclusion from ∼ 250 GeV until falling cross sections lead to a loss of sensitivity at ∼ 2.5 TeV, beyond which points are allowed. This exclusion limit falls with increasing tan β as the H + tb coupling is proportional to cot 2 β. 5 At moderate tan β the low mass exclusion limit shifts to H 0 /A 0 → 4b [180]. The kink at ∼ 100 GeV in the moderate tan β region is a result of B(H + → τ + ν τ ) falling to zero as B(H + → tb) rises to unity in this mass region, with B(H + → tb) ≈ 1 from m H ± m t . Sensitivity is lost from H 0 → as κ H 0 ∝ cot β, see Table 2. The flat cut off at 80 GeV results from H + → qq/τ + ν τ [176].
For the case of the 2HDM-II at low tan β, H 0 → γγ and H + → tb compete to give the most stringent exclusion bounds [179,181,182], with H + → tb being more sensitive at higher masses before the production cross sections fall at high masses. At moderate values of tan β, H + → τ + ν τ excludes masses up to ∼ 90 GeV, before the branching ratio for this channel falls, which leads to the allowed region at tan β of order 1 and masses of ∼ 100 GeV. This region is then ended by the H 0 → τ + τ − search in Ref. [183]. This search also gives the exclusion up to high masses in the large tan β region, exceeding the limits from H + → tb.
We are able to find lower mass bounds on the new Higgs bosons of 82 GeV and 86 GeV in the Type I and II 2HDM respectively, which are less stringent bounds than those found from the combined fit to flavour observables in Section 5.2. These bounds improve at low tan β in both models, and also at high tan β in the 2HDM-II, owing to the dependence of the couplings of these models to tan β. This is in line with the patterns we see in the flavour sector, where these couplings are also crucial.

Collider Outlook
Looking to the future, the LHC is due to be upgraded considerably, increasing the integrated luminosity to L HL-LHC = 3 ab −1 . This will have a sizeable impact, making it significantly harder for new particles to remain hidden. Here, we extrapolate the LHC data currently in HiggsBounds to this new luminosity by scaling the limits by a factor (L 0 /L HL-LHC ) for a search with a reference luminosity L 0 . There are important caveats to this extrapolation. One such caveat is that a number of searches are designed to precisely measure the SM Higgs behaviour. In these cases we match the bounds to what would be expected of the SM Higgs, as here we have a 125 GeV scalar in the exact alignment limit, which exactly matches the phenomenology of the SM Higgs. We make use of the SM Higgs predictions of Ref. [184] for this. Additionally, a number of searches included in HiggsBounds are from LHC run 1, during which the LHC operated at a centre of mass energy of √ s = 7 − 8 TeV, whilst the HL-LHC will operate at √ s = 13 − 14 TeV. We therefore need to reflect increased cross sections for these searches together with luminosity improvements (see also Ref. [185]). We again make use of the results from Ref. [184] to perform this extrapolation for SM Higgs searches. For BSM Higgs searches we use MadGraph5 aMC@NLO [175] to calculate the increase in the production cross sections at the higher centre of mass energy as a function of the BSM Higgs mass and scale the search data limits accordingly. We again perform a scan consisting of 50000 random points, the results of which are shown in Fig. 13, which includes the outline of the scan in Fig. 12 for ease of comparison. The permitted parameter space for the 2HDM is reduced at the future collider, based on extrapolating existing searches. In the 2HDM-I case this effect is only apparent at low tan β, as the sensitivity to new Higgs bosons is minimal at high tan β because the relevant couplings are proportional to cot 2 β. The lower limit in this region is from a LEP search [176], which we do not extrapolate here, hence the same cut off is present here as in Fig. 12.
There is no such issue in the 2HDM-II, and we find that the bounds improve significantly at the HL-LHC, bringing the lower mass bound into competition with those found from flavour constraints in Section 5.4, save for a small region with masses of ≈ 95 GeV and tan β ≈ 2.

Comparison with Flavour
In order to compare the bounds found above from collider searches to those we find from the flavour sector, we overlap results from both sections. No statistical combination between the two is attempted here as we do not use HiggsBounds to give a combined exclusion from all search data but only check if a point is excluded by any one individual search. We again present the results from the scenario in which all new Higgs bosons have degenerate masses and the alignment limit is exactly realised, meaning the combination presented in global fit in Ref. [13] for flavour results, where we simply extend these contours to lower charged Higgs masses to be more compatible with the collider results. We find that there is some degree of complementarity between the two sectors for the types of 2HDM that are examined here. In the 2HDM-I the LEP searches [176] that set a lower mass bound on the new Higgses outperform the exclusion from flavour observables, which lack sensitivity in the high tan β region. Conversely, in both cases, the flavour constraints following from precise measurements are more successful constraining tan β at at high masses, reflecting the loss of sensitivity in direct collider searches once the new particles are sufficiently heavy. Looking at the contours in Fig. 13, we can see that the extrapolation to the HL-LHC phase, which improves the sensitivity to new physics, enhances the complementary nature of the results in Fig. 14. In the case of the 2HDM-II the lower mass bound may exceed that determined from the flavour sector, depending on how future measurements in that sector line up with the scenarios outlined in Section 5.4.

A Note on Cosmological Context
As motivated by previous investigations [13,15,[186][187][188][189][190][191][192], we consider the possibility of generating a strong first order electroweak phase transition (SFOEWPT) through the 2HDM-I in order to fulfil the criterion for EW baryogenesis. Scanning large parameter regions, these find possibilities for a SFOEWPT to be achieved in the 2HDM-I, however this is when considering lower masses ( 1 TeV) than we find are favoured from our fits.
To evaluate the EWPT in the 2HDM-I at chosen benchmark points, we use the BSMPT package [193,194], where we refer to the package's documentation 6 for further information on the package, and we make use of it as discussed in Section 7.3 of Ref. [13]. In Table 5 we present the results for the strength of the EWPT at our selected parameter points. We first choose a benchmark of M = 50 TeV as a test of extreme masses, and then consider points motivated by our global fits in Table 3 as well as some variation in tan β. As the high masses of our best fits are suggested to be above the limit to achieve a SFOEWPT, we also consider a range of points with masses below 1 TeV, including the best fit point for the b → s observables. The alignment limit, cos(β − α) = 0, is taken for all benchmark points.
The results from these benchmark points are similar to those of the 2HDM-II. We find that the high mass scales around or above our best fit point (also enforcing close mass degeneracy) limit the strength of the phase transition. Sufficiently below our best fit points, where theoretical considerations now allow sampling with mass differences similar to those favoured in Ref. [188], we now find scenarios which support a SFOEWPT. In contrast to the 2HDM-II, these points with ξ c ≥ 1 are each within 1σ confidence of our best fit and are allowed by direct search data. In these scenarios, the global fit performs only slightly worse than the best fit points described in Table 3 with improvement over the SM at ∼ 1.5σ.

Conclusions
The search for new physics beyond the Standard Model is key to improving our understanding of physics at the smallest distances. Searches at the LHC for new interactions and states have so far not revealed any significant deviation from the SM expectation, while the precision of flavour and electroweak precision measurements corner the SM from different directions with some potential hints for LFU violation.
In this work we have performed a detailed investigation and comparison of the flavour and LHC exotics measurement programme, alongside its future extrapolation. We find that both flavour and current exotics searches can be significantly improved with increasing  Table 5. Table of results for the EWPT in the 2HDM-I. ξ c = ω c /T c is the parameterisation of the strength of the EWPT, with ω c the high-temperature VEV and T c the critical temperature. A SFOEWPT is indicated by ξ c > 1. In the limits chosen in Ref. [13] (and similarly discussed in Ref. [195]), λ 1 = λ 2 = m 2 h 0 /v 2 = 0.26 for all benchmark points, and so are not explicitly shown in the table. Additionally, λ 3,4,5 are not independent, but λ 3 + λ 4 + λ 5 = m 2 h 0 /v 2 = 0. 26. luminosity, thus allowing the experimental collaborations to cover a significant proportion of parameter space that is currently unconstrained. Extending the results of Ref. [13] for 2HDM-II constraints from indirect searches using electroweak precision, Higgs signal strengths, and flavour observables to the type I 2HDM, we find that the 2HDM-I statistically outperforms both the 2HDM-II and the SM in fits to the data. The best fit point for the 2HDM-I lies around The parameters of the 2HDM-I are found however to be much less constrained than in the 2HDM-II. An upper bound on the charged Higgs mass is only found at 1σ of 83 TeV, where for masses at or above ∼ 1 TeV, the additional Higgses are expectedly to be closely degenerate. For a lower charged Higgs mass, tan β is strongly constrained from below to be O(1) or higher, however this constraint is lessened with increasing charged Higgs mass. The alignment limit, cos(β − α) = 0, is favoured by our fits, however larger deviations than in the 2HDM-II are allowed, up to | cos(β − α)| = 0.4 at 5σ.
Expecting an increase in precision of flavour measurements from future colliders, we consider our flavour fits in several future measurements' projections, where these predictive fits encourage consideration of other models instead of (or in addition to) a 2HDM of type I/II to improve upon the SM. Furthermore, we consider the implications of these flavour fits on the electroweak phase transition of the early universe. It is found that the high masses encouraged by our fits do not support a strong EWPT for baryogenesis, however due to the 2HDM-I being less constrained than the 2HDM-II, sufficently low masses are possible within 1σ of our best fit.
Using data from direct searches for new fundamental particles at colliders we have been able to exclude large regions of the parameter space for the 2HDM of types I and II, highlighting possible areas to search for possible new physics. Looking forward, we have extrapolated a significant amount of data to match the expected performance of the HL-LHC and find that this further improves the constraints on new physics signals, to a level that can be competitive with indirect results from the flavour sector.
While some complementarity between flavour, Higgs physics and direct exotics searches remains, the finite energy coverage of exotics searches typically means a loss of direct LHC sensitivity for large masses approaching the decoupling limit.