Heavy Higgs Searches: Flavour Matters

We point out that the stringent lower bounds on the masses of additional electrically neutral and charged Higgs bosons crucially depend on the flavour structure of their Yukawa interactions. We show that these bounds can easily be evaded by the introduction of flavour-changing neutral currents in the Higgs sector. As an illustration, we study the phenomenology of a two Higgs doublet model with a Yukawa texture singling out the third family of quarks and leptons. We combine constraints from low-energy flavour physics measurements, LHC measurements of the 125 GeV Higgs boson rates, and LHC searches for new heavy Higgs bosons. We propose novel LHC searches that could be performed in the coming years to unravel the existence of these new Higgs bosons.


Introduction
Progress in physics often comes through understanding whether an intriguing observation is a relevant question worthy of further exploration. Particle physics is at a unique moment of its history and more than ever there is a need to identify the right questions: on one hand the Standard Model (SM) beautifully describes with an astonishing accuracy the vast majority of the data collected at various colliders; on the other hand, there are indisputable experimental evidences that this model has to be amended to account for instance of dark matter, to explain the matter-antimatter asymmetry, and also to be able to describe gravitational phenomena in the quantum regime. In addition, the values of the parameters defining the SM raise a number of questions: why is the Higgs boson mass so small compared to the scale of quantum gravity while it is subject to large quantum corrections? why is the electron mass 200 times smaller than the muon mass and 3,500 times smaller than the tau mass, while no quantum number distinguishes electrons from muons and taus?
The first question either requires new particles, new forces or new space-time structure, or large accidental cancelations in a multiverse. The second question is a less obvious guide toward the understanding of what should come next. Indeed, in the SM, an electron is different from a muon simply because it has a different Yukawa coupling to the Higgs boson. These Yukawa couplings are fundamental parameters of the SM and any values of these couplings are perfectly legitimate and not subject to large quantum corrections due to the chiral symmetry. But is it the correct description of nature? While the SM picture is compatible with the wealth of flavour measurements, the fact that the masses of the quarks and leptons truly originate from minimal Yukawa interactions with a minimal Higgs sector is still to be demonstrated. In particular, one needs to understand if the Yukawa pattern is completely arbitrary or if it is governed by some yet hidden structure(s). Additionally, the Yukawa couplings to first and second generations still need to be measured.
Within the minimal set-up of the SM, two particular empirical structures emerge: (i) the couplings of fermions to the Higgs boson are proportional to the fermion masses and (ii) these couplings are flavour diagonal. Were there new vector-like fermions mixing with the SM fermions, or were there more than one Higgs vacuum energies, these two emerging structures would be distorted and the Higgs interactions could trigger sizable flavour-violating processes. That is why the ∼2σ excess in the decay h → τ µ observed by the CMS Collaboration using LHC Run 1 data [1] raised much interest, although it has not been confirmed by subsequent analyses [2][3][4]. That is also why the search for t → ch decays is of high priority [5][6][7][8][9]. The confirmation of an excess or a null result would give invaluable information on the (non-)minimality of the Higgs sector and on the origin of flavour.
Conversely, the existence of sizable flavour-violating interactions in the Higgs interactions can have large impact on the phenomenology of additional Higgs bosons. This is the central question we want to explore in this paper. In particular, we demonstrate, by a concrete example, how the standard lower bounds on the masses of additional neutral and charged Higgs bosons can be easily evaded by departing from Type I and II two Higgs doublet models (2HDM) and by allowing for a small amount of flavour-violating interactions between the 125 GeV Higgs boson and the SM fermions at a level perfectly compatible with all flavour constraints. The flavour-violating window opens new rich phenomenological opportunities to which traditional searches for additional Higgs bosons remain blind. A full exploration of this new flavour territory is beyond the scope of this paper. Instead, we will limit ourselves to illustrating our points with a study of a 2HDM that does not fulfill the Glashow-Weinberg criteria [10] and can thus accommodate Higgs-mediated tree-level flavour-changing interactions. The most general form of this so-called Type III 2HDM has still a very large number of free parameters. 1 We will appeal to a specific texture in which the right-handed top quark is a singular player such that interesting signals can be expected in the forthcoming searches involving top quarks. The particular texture we will resort to was reconsidered recently by Ref. [12] to address the strong CP problem with an axion model free of the usual domain wall problem, making an interesting connection amongst three different topics (strong CP violation, dark matter abundance and flavour physics). Other scenarios have brought forward Higgs flavour-violating interactions to address some other experimental anomalies, notably the persistently troublesome value of the muon magnetic moment [13,14] or B → D ( * ) τ ν [15][16][17][18], B → K ( * ) µµ [19], or to explain the large mass hierarchy between the three generation quarks and leptons [20][21][22]. We leave these considerations aside, and focus our attention on how the presence of flavour-violating interactions impact the search for heavy Higgs bosons, and how new search strategies have to be designed to seize this novel opportunity.
Our paper is organized as follows. In Section 2, we present our extended Higgs model. In Section 3, we review the various experimental constraints from Higgs measurements and heavy meson decays. In Section 4, we advocate new collider signatures not yet explored at the LHC, as, for example, the associated production of a neutral or charged heavy Higgs boson followed by a subsequent flavour-violating decay, H → tc or H ± → cb. Finally, we present our conclusions in Section 5. Some useful formulae and details of the constraints on the model are collected in the appendices.

A top-philic flavour-violating two Higgs doublet model
A Type III 2HDM has generalized flavour-diagonal and off-diagonal Higgs couplings. In its most general form, the model is not consistent with experimental data on flavourchanging neutral current (FCNC) processes, since it predicts Higgs-mediated tree-level flavour transitions. However, as we will demonstrate, there exist specific flavour structures not based on Z 2 symmetries that are only weakly constrained by existing measurements (see also Refs.  for additional interesting models).
The most general renormalizable Yukawa sector involving two Higgs doublets, H u and H d respectively of hypercharge 1/2 and -1/2, and respecting the full local symmetry of the SM can be written as (withH i = iσ 2 Traditional Type II models correspond to vanishing non-holomorphic couplings,ˆ u,d, = 0, while Type I models have Y d, =ˆ u = 0. It is straightforward to write these Yukawa interactions in the mass-eigenstate basis of the Higgs bosons: 2) with v ∼ 174 GeV. Here α is the angle that defines the rotation from the interaction basis to the mass eigenstate basis of the CP even neutral Higgs boson fields. The angle β determines the rotation from the Higgs basis [42] to the interaction basis. In full generality, the fermion masses are diagonalized by bi-unitary transformations: In the mass eigenstate basis, the interaction Lagragian between the fermions and the Higgs bosons becomes It is well known [43] that in the decoupling limit, cos(β − α) → 0, the flavour-and CP-violating coupling of the light Higgs h, vanish but not those of H and A.
Tree level flavour violating couplings of the first two light generations of up type quarks in the Yukawa sector are tightly constrained by measurements of flavour violating processes. However, interesting phenomenology can be generated by allowing for flavour violating Yukawa couplings of the second and the third generations of fermions in the up sector in the light of recent experimental developments. For this we introduce a mixing between the right-handed charm and top quarks which can possibly be maximal since it involves only the right handed fermions. A similar argument can be placed for the lepton sector to motivate for possible interesting observations of flavour violating dynamics while suppressing flavour violation involving leptons from the first two families. The right handed mixing matrix for the up quark sector and the lepton sector are then given by For concreteness, we will further consider the set-up proposed by Chiang et al. in Ref. [12] that assigns a particular role to the top quark. This follows from a variant axion model endowed with a Peccei-Quinn symmetry acting only on the right-handed top. In this set-up, it follows that Upon performing the rotation defined in Eq. (2.4), we obtain the non-holomorphic couplings of the up-quarks and of the leptons: For the down quark sector, most flavour bounds can be respected by limiting the Yukawa couplings of Type II and hence we set d = 0. The explicit forms of the Higgsfermion couplings are listed in Appendix B. It is useful to note that the flavour-violating light Higgs couplings are proportional to the combination a = (tan β + cot β) cos(β − α) , (2.9) vanishing, therefore, in the cos(β − α) → 0 limit. For a matter of simplicity, in our phenomenological analysis of Sections 3 and 4 we will assume ρ u = ρ ≡ ρ. Relaxing this condition would lead to flavour-violating effects in the lepton sector which would be, to a certain extent, independent of the quark sector. Most notable, h → τ µ decays would be governed by ρ l while t → ch decays would be governed by ρ u rendering them theoretically uncorrelated. In addition, the scaling of the hτ τ coupling could, effectively, be different from the htt coupling in our model. While this could be an interesting study, we will not focus on it in this work. However, we will briefly comment on what changes when allowing different mixing angles in the up-quark and lepton sectors.

Constraints on the model parameters
In the model, due to the mixing with the second Higgs doublet, the couplings of the 125 GeV Higgs boson will be modified from what is predicted in the SM. We first analyze the constraints from Higgs coupling measurements in both flavour-conserving and flavourviolating processes (Sections 3.1 and 3.2, respectively). These bring about constrains on β − α, tan β, and ρ. In addition, the model allows for significant flavour violation, and therefore, its parameter space can be constrained by the measurement of several lowenergy flavour transitions, such as b → sγ, B → τ ν, R D , and R D * . We study the most significant constraints from low-energy flavour-violating processes on the parameter space of the model (Section 3.3). These processes add to the constraints on β − α, tan β, and ρ, and also constraint the mass of the charged Higgs (m H ± ), and, to a less extent, the mass of the neutral heavy Higgs bosons (m H , m A ). We comment on additional weaker flavour constraints in Section 3.4. Finally, we combine all relevant constraints in Section 3.5 and we demonstrate that heavy Higgs bosons as light as 200 GeV can be compatible with all constraints.
For combining the different constraints we use HEPfit [44], a code for the combination of indirect and direct constraints on High Energy Physics models. We use a Bayesian framework based on a Markov Chain Monte Carlo routine implemented within HEPfit. A custom version of the code is used since the code also allows for user-defined models. 2 We vary the parameters in the range: 0.1 ≤ tan β ≤ 15, −π ≤ ρ ≤ π, and 0 (β − α) π which are assigned flat priors. We restrict our investigation to not too large values of tan β, in anticipation of avoiding constraints from low-energy flavour-violating processes. The fourth parameter of concern, the charged Higgs boson mass m H ± , is varied in the range 200 GeV to 1200 GeV.

Higgs couplings measurements
The measurement of the effective couplings of the Higgs boson to SM particles puts strong constraints on the parameter space of the model, which predicts significant deviations of these couplings from the corresponding SM values at sizable values of the mixing angle, ρ, and cos(β − α) = 0. To study such constraints we use the results from the combination of measurements by both ATLAS and CMS collaborations in Run 1 [45]. Although more recent Run 2 ATLAS and CMS analyses are available [46-61], we decide to use the Run 1 combination, since the fits to and correlations amongst the ratios of the effective couplings are available, and this provides a more coherent way of comparing model parameters with experimental results.
In Table 1 we summarize the fit values and uncertainties for the effective couplings under the so-called "κ − λ framework" (coupling modifier ratio parameterization), along with the corresponding correlation matrix. These values have been obtained under the assumption that there is no exotic or invisible decay of the 125 GeV Higgs boson, in accordance with what is predicted by the model under consideration. Note that the asymmetric errors have been symmetrized since we use a symmetric Gaussian multivariate distribution to quantify the likelihood constructed from the fit to the experimental measurements. The

Constraints from the 125 GeV Higgs flavour-violating couplings
Measurements of the top-quark flavour-violating decays have been performed by both the ATLAS and CMS collaborations [7][8][9]. While for the branching ratio measurements of h → τ µ both CMS [1,4] and ATLAS [2] provide central values and RMS allowing us to build a likelihood profile, for t → ch we use only the ATLAS results [7,9] to put constraints on the model under consideration since the results from CMS [8] are cited only as upper bounds and have lower sensitivity than the latest ATLAS result. The numbers that are used in the fit are collected in Table 2, where we show a naive average of the BR(h → τ µ) and BR(t → ch) measurements. We plot these averages in Fig. 1 (dotted lines) in the BR(h → τ µ)−BR(t → ch) plane. 4 The solid lines represent the boundaries of the 68%, 95%, and 99% regions (green, yellow, and red line, respectively) of the prediction for BR(h → τ µ) and BR(t → ch) coming from the fit to the effective couplings in Table 1 only. The shaded region in Fig. 1 depicts how the measurements of the two flavour-violating processes constrain this prediction. From Fig. 1 it is clear that the experimental measurements of BR(h → τ µ) and BR(t → ch) constrain the parameter space of the model in addition to what is already constrained by 3 We have neglected the charged Higgs contribution to the gluon fusion production and the radiative decays, such that the Higgs coupling measurements directly constrain β − α, tan β and ρ independently of the value of the heavy Higgs masses. The charged Higgs contributions are, in fact, negligible in the entire parameter space. 4 The details of the analysis leading to this fit are reported in Appendix C. the Higgs effective couplings data. The combination of the flavour-preserving and flavourviolating Higgs data is presented in Appendix C. The parameter space of our model is further constrained by low-energy flavour observables, as discussed in the next section.

Low energy flavour constraints
Having flavour off-diagonal couplings, the several Higgs states can leave very distinctive signatures in FCNC and tree-level charged-current processes in low-energy mesonic decays.
To validate the parameter space of our interest, it is important to look at the possible constraints coming from these decays both in inclusive and exclusive channels, as well as from neutral meson oscillations. As we will discuss, flavour constraints favor the part of the parameter space that involves a relatively low value of tan β. While the Higgs couplings data and the flavour-violating processes discussed in Sections 3.1 and 3.2 put bounds on only tan β, ρ and β − α, the low-energy flavour-violating processes that we consider here also put constraints on the mass of the charged Higgs boson. As we will comment in Section 3.4, constraints on the mass of the heavy neutral Higgs boson are much milder. In a Type II 2HDM, the bound on m H ± coming from the measurement of BR(b → sγ) is very stringent [62,63] and is in the 570-800 GeV range at 95% confidence level 5 , independent of the value of tan β, for tan β 2. However, this constraint can get weakened in a Type III 2HDM [16] because of a possible destructive interference with the SM contribution for certain ranges of u 23 and u 32 allowing for much lower values of m H ± to remain compatible with measurements. In terms of the parameters in this model, a non-zero ρ can bring about contributions that interfere destructively with the SM contributions allowing for compatibility with experimental measurements for relatively low values of m H ± . The tan β dependence of the new physics (NP) contributions from our model is also quite different from a Type II 2HDM, since both u i L d j R and u i R d j L charged Higgs couplings are tan β-enhanced for ρ = 0 (see Eqs. (A.3)-(A.4)).
In addition to b → sγ, the charged-current process B → τ ν can get contributions due to a flavour-violating charged current mediated by the charged Higgs boson. Finally, the  Table 3. The experimental measurements for the flavour changing observables used in the fit. The statistical and systematic uncertainties on the R D and R D * measurements have been summed in quadrature. These two measurement have a correlation coefficient of −0.23.
observables R D and R D * involving a b → cτ ν transition can also set constraints on our parameter space. In our analysis, we do not try to explain the long-standing discrepancy between the SM predictions and the experimental measurements of these observables, but we simply include the two observables in a global fit. We collect the measurements and SM predictions of the low-energy flavour-violating processes that we use in this study in Table 3. More details about the decay modes b → sγ, B → τ ν and the observables R D and R D * can be found in Appendix D where we also discuss the subtleties involved in performing the necessary calculation in the particular model under consideration.
We note that additional flavour observables will be modified by the charged Higgs exchange. Examples are leptonic meson decays, However, the constraints coming from the measurement of these transitions are always much milder than the ones arising from b → sγ and B → τ ν. For this reason, we do not introduce these additional observables in our global fit.

Other indirect constraints
Given the flavour structure of the model, one could expect that the amount of flavour violation encoded by the parameter ρ can be severely constrained by additional flavour observables beyond those that we have discussed in Section 3.3. While we do not present a detailed analysis of additional constraints, in this section we discuss why FCNC observables do not affect our current study. In particular, we will discuss why neutral Higgs mediated processes are only weakly constrained, leading to very mild bounds on the heavy Higgs mass, m H , in the regions of tan β, ρ allowed by the charged current constraints discussed in the previous section.
In generic Type III 2HDMs, amongst the tree-level FCNC processes mediated by the several neutral Higgs states, important constraints can arise from ∆F = 1 leptonic meson Neutral Higgs mediated tree-level contributions to lepton flavour-violating (LFV) decays like τ − → µ − µ + µ − and τ − → e − µ + µ − arise in our model. However, the contributions are rather small, since they are proportional to the lepton mass of the respective generation, and the experimental bounds are relatively weak. This translates into weak bounds on 23,32 , 13,31 which do not affect our analysis. Better measured LFV processes like µ − → e − e + e − do not receive tree-level NP contributions since 12,21 = 0 in our model. The bounds from radiative LFV processes like τ − → µ − γ and τ − → e − γ are significantly weaker.
Additional low energy observables will be affected by the exchange of neutral Higgs bosons. There is a long-standing discrepancy between the SM prediction and the measurement of the muon magnetic moment, a µ = (g − 2)/2. In our model, loop contributions mediated by the neutral (and charged) Higgs boson exchange are generated. As shown by the couplings in Eqs. (A.5)-(A.6), a very large value of 22 (tan β + cot β) is needed to obtain a relatively sizable NP effect. Since 22 ∝ m µ /v, a µ receives only a small NP effect in the region at low values of tan β, as required by the measurement of b → sγ.
Finally, also electroweak precision observables like the S, T oblique parameters can set important constraints on the parameter space of a 2HDM. These electroweak parameters were studied in a generalized Type III 2HDM e.g. in Refs. [64,65]. In these studies, it has been shown that the T parameter is generically the most important constraint, and its bound depends on the mass splitting of the charged Higgs boson and the most massive neutral Higgs boson. Since, in our investigation we do not fix the relative mass of the heavy Higgs states and they can possibly be quasi-degenerate, the T parameter will not put significant constraints on the parameter space under consideration.

Combining constraints
As a final step to determine possible benchmark points that are allowed by all the constraints, we perform a combined fit using the Higgs effective coupling results tabulated in Table 1, the Higgs and top flavour-violating decay results collected in Table 2, and the charged-current process measurements listed in Table 3.
We show the results of our combination in Fig. 2, as a function of the free parameters, β − α, ρ, tan β, m H ± . As we have deduced before, the primary constraints come from measurements of Higgs coupling, of the branching ratios of t → ch, and of b → sγ. The favored values of β − α is clustered at around π/2, corresponding to the decoupling or alignment limit. The dashed contours in Fig. 2 represent the 68%, 95%, and 99% contours from constraints coming only from the Higgs flavour-conserving and flavour-violating measurements discussed in Sections 3.1 and 3.2 (see also Fig. 10 in Appendix C). As can be seen from the right panels of Fig. 2, the addition of the flavour violating low-energy observables, and in particular of b → sγ, brings about a preference for lower values of  Table 1), h → τ µ, t → ch (see Table 2) and flavour measurements (see Table 3). The green, yellow, and red regions are the 68%, 95%, and 99% regions respectively. The dashed contours represent the 68%, 95%, and 99% contours from constraints coming only from the Higgs flavour-conserving and flavour-violating measurements. The black dot in each plot marks the benchmark point we will use for the discussion of the collider phenomenology in Sec. 4, and they correspond to ρ = 1, tan β = 1, and β − α = 1.4454 (cos(β − α) = 0.125). tan β as well as a bit smaller values of |ρ|. The dramatic change of the constraints to the parameter space can be seen comparing dashed and solid lines in the lower right panel of To complete this analysis, we choose a benchmark point that we will use in the next section to discuss the collider phenomenology of our model. The point of our choice is marked with a black dot in the panels in Fig. 2, and it is in the 68% probability region in all the 2D marginalized posterior distributions and corresponds to ρ = 1, tan β = 1, Figure 3. Some of the most important Feynman diagrams contributing to heavy Higgs production in our model (first two panels), as compared to a Type II 2HDM (last two panels). We do not report the gluon fusion diagram since it is common to the several types of 2HDMs.
In what follows we shall see that this point in the parameter space brings about some interesting phenomenological implications for heavy Higgs boson production and decay. As we will show, varying the point in the favored region of parameter space will not qualitatively affect the heavy Higgs boson phenomenology.

Collider signatures
One of the primary goals of our study is to shed light on possible interesting collider signatures that have not been searched for yet at the LHC. As we will show in this section, over a broad region of the allowed parameter space of the model, the phenomenology of the heavy Higgs bosons, H, A, H ± , differs significantly from what is typically assumed in direct searches at the LHC. Particularly, the heavy Higgs bosons can have large branching ratios for flavour-violating decays, similarly to what is predicted by different flavour structures such as e.g. the flavourful 2HDM (F2HDM) of Refs. [20][21][22]. Furthermore, also some of the main production mechanisms will be different than the ones predicted by a Type II 2HDM (see Fig. 3 for the most important Feynman diagrams contributing to the heavy Higgs production in our model, as compared to a Type II 2HDM. In the figure, we do not report the gluon fusion diagram since it is common to the several types of 2HDMs). For definiteness, and following the discussion in Section 3, we consider a benchmark scenario with tan β = 1 and cos(β − α) = 0.125, and study how the phenomenology depends on the assumed value of ρ and the heavy Higgs boson masses. The choice of a different benchmark with cos(β − α) = 0 would not alter the phenomenology qualitatively.
The predicted cross sections and branching ratios were obtained using FeynRules v2.3 [66,67] and MadGraph 5 v2.3.3 [68]. In order to obtain next-to-leading-order (NLO) cross sections from MadGraph, we implemented our own model in FeynRules based on the NLO implementation of the 2HDM [69] and modifying it to a 3, 4 or 5 massless flavour theory as and when necessary. The default MadGraph run cards were used with the NNPDF2.3 PDF sets [70] with matched order. To study the gluon fusion production cross sections of the neutral Higgs boson we used HIGLU [71], since this code allows for the independent extraction of the top-quark and b-quark loop contributions and their interference. We then rescaled these cross sections with the modified couplings from our model.

Phenomenology of the neutral heavy Higgs boson
As shown in Section 2, the couplings of the neutral heavy Higgs bosons to fermions depend strongly on the value of ρ. In the case of ρ = 0, the couplings are purely flavour conserving, whereas for ρ = π/2 the flavour violating couplings can be even larger than the corresponding flavour-conserving ones.
To illustrate the effect of a sizable value of ρ (which is nevertheless consistent with the experimental constraints from Higgs and flavour transition measurements), in Fig. 4 we present the branching ratios (left panel) and production cross sections at the 13 TeV LHC (right panel) for the CP-even neutral Higgs boson as a function of its mass, having fixed ρ = 1 (solid lines in the plots). 6 For comparison, we also show the corresponding curves in the flavour conserving case ρ = 0 (dashed lines in the figure). In this limit, the Higgs couplings to third generation show a Type IV structure, namely Type II-like for the bottom and top couplings and Type I-like for the tau coupling.
As we will also see in Fig. 5, sizable values of ρ result in H → tc being the main decay mode. The branching ratio into tt is suppressed by about a factor of ten, when compared to the well-studied Type II 2HDM at low tan β (see the dashed yellow line in the figure). The choice of tan β = 1 leads to only very small branching ratios into bb and τ + τ − , as it can be seen from the expression of the couplings in Appendix B. Particularly, the τ + τ − decay mode is quite suppressed if compared to the Type IV 2HDM prediction (see solid vs. dashed red lines in the figure). As our benchmark point requires cos(β − α) = 0.125 and thus requires some mixing between the two Higgs doublets, the branching ratios into W + W − and ZZ are non-zero but small, in the range of ∼0.5%-10%, depending on mass. Finally, the branching ratio for the flavour-violating decay H → τ µ is generically larger than the one for H → τ τ , even if it is, in any case, quite small (see e.g. Refs. [73][74][75] for models predicting sizable lepton-flavour-violating branching ratios).
The main production mechanisms are gluon fusion, followed by associated production with a top quark and a c-quark (tcH). This novel production mode has a sizable cross section due to the sizable Htc coupling (see Eq. (A.11)). The productions in association with top quarks (ttH) and b-quarks (bbH) are relatively small, both having cross sections of the order of tens of fb for a heavy Higgs boson with a mass of a few hundred GeV. Particularly, the ttH cross section is rather suppressed if compared to a Type II 2HDM (solid vs. dashed yellow lines in the figure); the bbH cross section is the same as in a Type II 2HDM, since the down quark sector of our framework has the same coupling structure as a Type II 2HDM.
In Fig. 5 [88,89] do not set relevant constraints on the parameter space, and are not shown in the figure. Figure 5 also shows the 68% and 95% favored regions that are obtained from combining all the constraints listed in Tables 1, 2 and 3, setting tan β = 1 and cos(β − α) = 0.125, and marginalizing over the value of m H ± (shaded green and yellow regions, respectively). The dashed lines correspond to the same regions but assuming The processes with the largest signal rate at sizable values of ρ (∼1) are pp → H(→ tc) and pp → tcH(→ tc), with cross sections that can reach ∼ 3 pb and ∼ 1 pb, respectively, in the region of parameter space not yet probed by the present LHC heavy Higgs searches (blue shaded region) and in agreement with Higgs and flavour measurements (green and yellow shaded regions). These two processes are also complementary in their coverage of ρ, with the former process peaking at |ρ| ∼ 0.6 while the latter peaks at |ρ| ∼ 1.5, for the values of the other parameters, tan β = 1 and cos(β −α) = 0.125, chosen in our benchmark. The additional production modes, bbH and ttH, lead to cross sections that are at least one order of magnitude smaller than tcH and pp → H, and, therefore, are not the first smoking gun of our framework.
Experimentally the pp → H(→ tc) process is quite challenging because of the relatively small signal rate in comparison with the large expected backgrounds, mainly coming from the tt+jets and W +jets SM processes. Also, for heavy Higgs bosons with masses of few hundred GeV, the resulting top quarks would typically not be boosted, which experimentally implies focusing on leptonic top-quark decays in order to be able to trigger on the events. 7 Here we consider only those searches that take into account the interference effect between the nonresonant tt background and the signal, which significantly distorts the spectrum, creating a peak-dip structure [81][82][83][84][85][86]. A similar signature is probed by searches for pp → H + (→ tb) [90] and pp → W (→ tb) [91-93], which have some difficulty in probing low values of mass (< 500 GeV), even with two b-tagged jets in the event, a key requirement to suppress backgrounds from W +jets production. In the case of our signal of interest, one of those bottom jets would be replaced by a charm jet, which is harder to discriminate from light-flavour jets, although the LHC experiments are making progress on optimizing the performance of charm-tagging algorithms [94]. It will be important to assess the prospects for probing the process pp → H(→ tc) in the coming years of the LHC using these new developments, even if we expect that reaching the required level of sensitivity would be rather challenging.
On the other hand, the process pp → tcH(→ tc), with two top quarks and two additional charm jets in the final state, provides more experimental handles to suppress the backgrounds. There are two experimental signatures that are potentially promising. (1) In the case of only one top quark decaying leptonically, the final state is characterized by one lepton, some missing transverse momentum (E miss T ) due to the presence of a neutrino, and at least six jets, of which two are bottom jets and two are charm jets. In this case the main background is tt+jets production, and in particular tt+bb and tt+cc, which are subject to large uncertainties both theoretical and experimental. However, the LHC experiments have developed sophisticated search strategies exploiting large-statistics subsidiary data samples for the purpose of constraining background uncertainties. Examples are searches for tbH ± (→ tb) [95] and bbH(→ tt) [96], although they are not optimized for this signal. A natural evolution of those searches to target this signal would include exploiting kinematic reconstruction of the heavy Higgs boson mass, and possibly using explicit charm tagging as a way to suppress the tt+bb background. It will be very interesting to assess the reach of this possible search. (2) An additional and even more promising avenue to probe the pp → tcH(→ tc) signal is through the final-state signature of same-charge dilepton plus bottom and charm jets, since half of the events have two same-charge top quarks. In this case SM backgrounds are heavily suppressed compared to the lepton-plusjets final state discussed previously, although instrumental backgrounds from tt+jets with charge-misreconstruction or non-prompt/fake leptons need to be precisely estimated using data-driven techniques. Examples of this kind of searches, currently not tailored for the pp → tcH(→ tc) signal, but which could be optimized accordingly, are searches for pair-production of vector-like quarks [97] or anomalous four-top-quark production, e.g. via pp → ttH(→ tt) [97, 98] (see also Refs. [82,84,99]).
Finally, searches for same-charge dilepton or trilepton signatures are also very promising to probe the pp → tcH(→ tt) process, with three top quarks in the final state, provided the H → tt branching ratio is sufficiently large. 8 However, we do not discuss this channel in great detail here since the flavour violation in the model under consideration leaves more dominant signatures in the decay modes discussed above.

Phenomenology of the charged Higgs boson
The phenomenology of the charged Higgs boson also changes dramatically at non-zero values of ρ, if compared to the better studied Type I-II 2HDMs. In Fig. 6 we present the branching ratios (left panel) and production cross sections (right panel) for the charged Higgs boson as a function of its mass, assuming cos(β − α) = 0.125, tan β = 1, and ρ = 1 (solid lines) or ρ = 0 (dashed lines). Note that, contrary to the neutral Higgs case, in the charged Higgs case the ρ → 0 limit correspond to formally finite but numerically small deviations from Type II Yukawa couplings in the quark sector. In the lepton sector, the couplings do not match any of the Type I-IV 2HDMs. For the ρ = 1 benchmark, the H ± → cb decay mode becomes dominant, with a branching ratio well above ∼70%, followed by H ± → tb with a branching ratio below ∼20% and quite smaller than the corresponding prediction in a Type II 2HDM (see solid vs. dashed red lines in the figure). This exact ratio of branching ratios is due to the specific benchmark chosen. Particularly the ratio of the tb and cb couplings can be approximated by ∼ (1+tan showing that the hierarchy between the tb and cb decay modes can also be reversed, even if typically the cb coupling is larger than the tb one. Because of the non-zero value of cos(β − α) the branching ratio H ± → W h is also sizable. Furthermore the very well studied τ ν mode is very suppressed and has branching ratios at the level of few ×10 −5 . This branching ratio is smaller than the corresponding one in a Type II 2HDM and in the ρ = 0 case (solid vs. dashed yellow lines in the figure). For our benchmark, because of the enhancement of the H ± cb coupling relative to the H ± tb coupling, the main production mechanism is pp → cbH ± , with a cross section orders of magnitude above that for pp → tbH ± , and much larger than what is predicted by a Type II 2HDM (solid vs. dashed blue lines in the figure). The H ± tb cross section is, instead, a factor of a few smaller than in a Type II 2HDM (solid vs. dashed red lines in the figure). As a result, in our scenario the usual signal process searched for at the LHC, pp → tbH ± (→ tb) [95, 101], is quite suppressed, while there are other signatures that are much more promising and remain relatively unexplored.
In Fig. 7 we show σ × BR for several interesting H ± at the 13 TeV LHC, as a function of ρ and of the charged Higgs boson mass, m H ± . The 68% and 95% favored regions that are obtained from combining all the constraints discussed in the previous sections and listed in Tables 1, 2 and 3 are shaded in green and yellow, respectively. Over most of the parameter space, except when ρ is small, the dominant process is pp → cbH ± (→ cb), with σ × BR ranging from O(500 pb) at m H ± ∼ 200 GeV to O(1 pb) at m H ± ∼ 800 GeV. This process can, in principle, be probed by inclusive searches for dijet resonances at the LHC. Currently a few searches probe resonance masses below 1 TeV (see Refs. searches at √ s = 13 TeV), which is the mass range of interest in this study. In particular, among these searches, the most stringent constraint on our model is set by Ref. [102], and is displayed as a shaded blue region in Fig. 7. To obtain this region, we recast the ATLAS analysis [102]. We use MadGraph in the four-flavour scheme to compute cross sections and generate parton-level events at √ s = 13 TeV for the process pp → cbH ± (→ cb), for different values of the charged Higgs mass. Subsequently, we require the two leading jets to have a p T larger than 185 GeV and 85 GeV, respectively, and |η| < 2.8, as required by the ATLAS analysis. The requirement for |y * | < 0.3 or |y * | < 0.6, depending on the signal region, is also added (y * is the difference in rapidity of the two leading jets, y * ≡ (y 1 − y 2 )/2).
Finally, we extract the bound on the parameter space, by comparing the cross section of our model after cuts, with the excluded cross section. As can be appreciated, most of the mass region is still allowed even at high values of ρ. However, these analyses only use a fraction of the Run 2 dataset and further optimizations for this signal should be possible, including making b-(and c-)tagging requirements. Therefore, we tentatively conclude that dedicated future searches for resonances in multijet final states would be quite important to probe a significant fraction of the parameter space in this scenario.
The process with the next highest rate is pp → cbH ± (→ tb), with a σ × BR of several pb below a mass of ∼600 GeV. Since there is a high probability that the associated bquark and c-quark fall outside the ATLAS and CMS detector acceptance (only a few percents of events have both charm and b-jets satisfying the standard p T > 20 GeV and |η| < 2.5 requirements) the main signature is a tb resonance, for which there are dedicated searches [90-93]. In practice, however, current searches [93] at √ s = 13 TeV are able to probe σ × BR ∼ 1 pb at a mass of 1 TeV, about one order of magnitude above the expected σ × BR in our scenario. With sufficient integrated luminosity and further optimisations, these searches may eventually be able to probe our scenario at masses lower than 1 TeV, although it is not clear if masses as low as 200 GeV can be reached.
Other signals have smaller cross sections, with σ × BR exceeding 0.5 pb only at low mass. This includes not only the standard LHC search mode, pp → tbH ± (→ tb), but also new processes such as pp → cbH ± (→ W ± h) and pp → tbH ± (→ cb). In the case of pp → tbH ± (→ cb), the final state signature is the same as in LHC searches for tt → W bH ± (→ cb)b [108], which so far have assumed a light charged Higgs boson in top-quark decays and, thus, are currently not optimized to probe this signal. Our model offers, therefore, a motivation to search for H ± → cb resonances even above the top threshold. 9 In the case of pp → cbH ± (→ W ± h), since the associated b-quark and c-quark are often not reconstructed, the signal is similar to the one probed in searches for the SM Higgs boson in the W h production mode, except for the presence of a W h resonance. Existing searches for W → W h(→ bb) resonances at √ s = 13 TeV [109,110] currently target resonance masses above at least 500 GeV, although masses down to 300 GeV were also probed by Run 1 searches [111]. The most significant constraint on the parameter space of our model is set by the search in Ref. [109], which leads to the bound shaded in gray in Fig. 7. Due to the large cbH ± production cross section, 10 already sizable regions of parameter space 9 See also Ref. [22] for another model predicting this signature. 10 Note that the cbH ± production cross section in our model can be much larger than the tbH ± production are probed, leading us to conclude that there will be very good prospects for W h searches to probe our model, in the case of cos(α − β) = 0. Modifications of the current searches can also be envisioned, particularly targeting the cbH ± associated production through the requirement of an associated heavy-flavour jet, which may lead to further improvements in sensitivity.
In the model used in this work the branching fraction of H ± → cs is negligibly small compared to the other decay modes in contrast with come other models where H ± has significant flavour violating couplings with the second generation of quarks. A search by CMS for this decay channel can be found in Ref. [112].

Comparison with other top-philic flavour-violating 2HDMs
As we have demonstrated in this section, 2HDM frameworks that single out the third generation of quarks and leptons have a very different phenomenology, if compared to 2HDMs with natural flavour conservation. Particularly, it is very easy to evade the constraints from the present LHC direct searches for new Higgs bosons and, at the same time, predict sizable cross sections for signatures not yet searched for at the LHC. Different scenarios will, however, predict a different pattern of branching ratios and cross sections. In this section, we discuss how to disentangle experimentally between the framework discussed in this paper and other scenarios, in particular the F2HDM proposed in Refs. [20,21] that generates the third generation fermion masses mainly through the SM Higgs mechanism and the 1st and 2nd generation masses through the Higgs mechanism of an additional Higgs doublet.
Both frameworks predict sizable branching ratios for H, A flavour-changing decays. However, in the F2HDM the flavour-violating couplings between 2nd and 3rd generations scale as m 2 nd /v tan β, where m 2 nd is the mass of the second generation fermions, and have therefore a similar scaling as the flavour-diagonal couplings to second generations. In our framework, instead, the 2 ↔ 3 flavour-violating couplings in the up and lepton sectors can be quite larger than the coupling to second generations, since c 23 ∼ (m 3 rd /v)(cot β + tan β) sin ρ, where m 3 rd is the mass of the third generation fermions (see Eq. (A.11)). Similarly, flavour-violating couplings involving the first generation quarks and leptons are suppressed by the first generation masses in the F2HDM; in the minimal realization of our framework where the third generation up quark and lepton mixes only with the second generation, there are no neutral Higgs couplings to first and second generations. If we consider a richer flavour structure that mixes the first generation quarks and leptons with the second and third (see the discussion in Section 2), the two flavour-violating couplings will scale as c 13 ∝ m 3 rd /v and c 12 ∝ m 2 nd /v, eventually generating larger flavourchanging effects that, however, will be constrained by low energy flavour observables like D-D mixing.
This different scaling of the several couplings results in a quite different phenomenology for the heavy Higgs bosons. In particular, searches for di-muon and di-jet resonances in a Type II 2HDM, even at very large values of tan β, and can lead to a much larger number of charged Higgs bosons produced at the LHC.
will be more relevant to test the neutral Higgs bosons arising in the F2HDM, than those of the framework discussed in this paper, where the only relevant bounds arise from searches for H → ZZ (see Fig. 5). Similarly, also the phenomenology of the charged Higgs boson presents both differences and similarities with the one predicted by the F2HDM. In both 2HDMs the flavour-violating decay H ± → bc has sizable branching ratios. However, contrary to the F2HDM, in our framework this decay mode is typically the dominant one together with the tb decay mode, and has a larger branching ratio than the cs decay mode. Additionally, the charged Higgs boson is typically more copiously produced at the LHC in our framework than in the F2HDM, due to its larger flavour-violating cb coupling, that scales as m t /v.
Another class of models that warrants a mention here is the one proposed by Branco-Grimus-Lavoura (BGL) in Refs. [25,113,114]. In these models the flavour-changing couplings of the neutral Higgs sector are connected directly with the quark and lepton mixing matrices. Recently, these models have been used to study signatures like h → τ µ and t → ch also considered in our work. It has been shown that these models can saturate the experimental bounds on these channels [30]. However, in these models the flavour offdiagonal couplings of the neutral Higgs are proportional to the products of the off-diagonal elements of the mass mixing matrices which, in the quark sector, are suppressed by powers of the Wolfenstein parameter. Hence, the flavour-violating effects can never be as large as those allowed in the model we use in this study, where we have shown that the parameter space of the model is actually significantly constrained by the experimental measurements of h → τ µ and t → ch.
Flavour violation in the top sector can also be enhanced while suppressing the one in the down sector by assuming the Cheng-Sher ansatz [23]. The Yukawa coupling structure is assumed to have the form , where i, j denote the family indices and λ ij can be of O (1). This makes the flavour violating effects proportional to the geometric mean of the masses and hence suppressing the FCNC effects in the first two families while allowing significant effects in the processes involving the third family, especially the top quark. The flavour violating effects can be significant in the top sector, a study of which can be found in Ref. [17]. However, since their study focuses on the alignment limit, where the SM Higgs sector is decoupled from the NP Higgs sector, they do not predict any NP effects in channels like t → ch or h → τ µ, in contrast with the model that we study in this work.

Conclusions
The discovery of the Higgs boson came with pride and prejudice. Pride that a new dynamical principle, namely spontaneous symmetry breaking, envisioned by theorists about 50 years ago, found its incarnation in Nature. Prejudice that the Higgs boson is not the last particle of the SM but rather the first particle en route towards physics beyond the Standard Model. The Higgs boson remains so far the only particle that exists in a single species. Of course, this could well be the result of large quantum corrections that push to the mass of elementary scalars to high value. If they are present around the weak scale, extra Higgs bosons are very much expected to be accompanied with a plethora of new particles, be they supersymmetric or composite, that would balance the large corrections to the Higgs masses. These new particles can affect the production and the decay of these extra Higgs bosons in a significant way and it is therefore difficult to put robust and model-independent bound on the heavy Higgs boson masses. Therefore the stringent bounds obtained in standard searches in minimal Type I and Type II 2HDMs should be taken with care.
The purpose of this paper has been to show that simple departures from minimal setups fool the standard direct and indirect constraints from Higgs coupling measurements and flavour data. Rich opportunities for discovery open up in unexplored channels. Not only neutral and charged Higgs bosons as light as 200 GeV can exist, but they can be copiously and predominantly produced and decay in channels that have not be explored yet. For instance, we showed that a neutral heavy Higgs boson can be produced in association with a top and charm quarks and later decay into another top-charm pair, leading to a final state with same charge dilepton and bottom and charm jets that can be easily emerged from a not so dominant background. A heavy charged Higgs boson can be produced in association with a bottom and charm quarks and can decay into another bottom-charm pair with a total rate at or above 100 pb. It will be very interesting to extend the present LHC program for searches of new Higgs bosons, to include the plethora of new signatures predicted by our model (see Table 4). For sure rather spectacular signatures are expected and wait for the interest of the experimental community to reveal the first direct evidence of new physics and to unravel the origin of flavour, which remains one of the deepest questions of high-energy physics.
de Economía y Competitividad under projects FPA2015-69260-C3-1-R and Centro de Excelencia Severo Ochoa SEV-2012-0234. A.P. would like to acknowledge support from the ERC Ideas Starting Grant n. 279972 "NPFlavour". S.G. is grateful to the hospitality of the Kavli Institute for Theoretical Physics in Santa Barbara, CA, supported in part by the National Science Foundation under Grant No. NSF PHY11-25915, the Aspen Center for Physics, supported by the National Science Foundation Grant No. PHY-1066293, and the Mainz Institute for Theoretical Physics (MITP) where some of the research reported in this work was carried out. Finally, the authors would like to thank the Galileo Galilei Institute in Florence, the Centro de Ciencias Pedro Pascual in Benasque, and the Mainz Institute for Theoretical Physics, where part of this project has been conducted during the workshops "Gearing up for LHC13", "Higgs Tasting" and "The TeV scale: a threshold to new physics?", respectively.

A Higgs couplings to fermions
The physical Higgs-quark and Higgs-lepton couplings (H 0 k = H, h, A) [115] are given by: 11 For the coupling of the charged Higgs boson with leptons, we sum over the neutrino flavour (ν L ≡ f ν f L ). We have defined the coefficients x k q as Using the expressions in (2.8) for u ij , we find that the 125 GeV Higgs couplings to up quarks are given by For completeness, we also report here the corresponding couplings of the heavy Higgs bosons, H and A: tan β (for f = u) . (A.10) The flavour-violating couplings of the several Higgs bosons with charm and top quarks are given by The couplings of the Higgs bosons to leptons have a similar form, requiring an exchange of ρ u → ρ and the τ (µ) coupling being similar to the top(charm) coupling. Finally, all couplings of the Higgs bosons to down-type quarks are as in a Type II 2HDM.

B Constraints from Higgs effective couplings
The Higgs effective couplings that parametrize the fit to the measurements of Higgs boson productions and decays can be written in terms of the reduced couplings, κ, as  Table 1. We vary the parameters in the range: 0.1 ≤ tan β ≤ 15, −π ≤ ρ ≤ π, and 0 (β − α) π, with flat priors. The green, yellow, and red regions are the 68%, 95%, and 99% favored regions, respectively.
In all generality, in the absence of exotic Higgs decays, the scaling of the light neutral Higgs boson width is given by [45]: where we are neglecting the decays to the light generations, s, u, d, e. In the case of no beyond the Standard model particles running in the loops for the Higgs to gluon, photon and Zγ effective couplings, the scalings of the couplings of the Higgs to the gauge bosons are given by: and the scalings of the fermionic couplings, κ f , can be derived using for f = t, b, τ and c h f couplings given in Appendix A. It should be noted that, in our work, we assume the same flavour structure in the up quark sector and in the lepton sector and, therefore, κ τ = κ t .
In Fig. 8 we present the 2D marginalized posterior distributions for the parameters ρ, tan β, and β − α, as constrained by the Higgs effective couplings data listed in Table 1. As expected, the value of |β − α| is constrained to be not too large, since the 125 GeV Higgs boson has SM-like properties. The values for ρ and tan β are, instead, relatively unconstrained by the fit. Our analysis leads to results that are qualitatively similar to those of Ref. [12], where the same model was used for the analysis. However, it should be noted that, in contrast to Ref. [12], we use the final results from the Run I Higgs coupling combination, taking into account correlations. Furthermore, our statistical approach is different, as we obtain the 2D plots by marginalizing over the third parameter rather than assuming a fixed value for it. 12

C Constraints from flavour-violating decays involving the Higgs boson
The off-diagonal fermion couplings of the light neutral Higgs boson are, in general, non-zero, and can lead to a measurable value for branching ratios of h → τ µ and t → ch. 13 Allowing for deviations in the production cross section of the light Higgs boson, the experimental measurement of the h → τ µ branching ratio should be compared to [12] where we denote with σ pp→h and BR th (h → τ µ) the production cross section of the 125 GeV Higgs and the branching ratio of its decay to τ µ, as predicted in the scenario we study. For simplicity, we have only considered the gluon fusion production mechanism and the decay into bottom quarks, W W , ZZ, gluons, and tau leptons, leading to the terms in the denominator. The other decay modes will only bring small corrections that can be safely neglected. In the rest of the text we will refer to BR exp (h → τ µ) as BR(h → τ µ). The quantity a is defined in Eq. (2.9).
In the quark sector, the branching ratio for the decay t → ch can also be sizable, and can be written as [12] 14 Other recent studies of the decay t → ch in a general 2HDM framework can be found in Refs. [30,[116][117][118] and those for h → τ µ can be found in Refs. [18,[118][119][120][121][122][123][124][125]. In Fig. 9, we show the individual constraints from BR(h → τ µ) and BR(t → ch) in the ρ vs. β − α plane, having marginalized over the other free parameter, tan β. The green, yellow, and red regions are the 68%, 95%, and 99% favored regions, respectively. The left panel corresponds to the bounds on the parameter space coming from t → ch data. The right panel of Fig. 9 shows the constraints coming from the h → τ µ decay. The combination of the two constraints is shown in the middle panel. As can be gauged from 12 The results in Ref. [12] are presented in terms of tan β, ρ and a (defined in Eq. (2.9)). 13 A richer flavour structure including mixing of the third generations, not only with the second one, but also with the first generation would lead to additional interesting decay modes such as h → τ e.
14 As shown by the structure of the Higgs couplings presented in Section 2, the decay t → uh does not occur due to the absence of the corresponding flavour-violating couplings.   Table 1) and the flavour-violating h → τ µ and t → ch decays (see Table 2). The dashed lines show the constraints from only the Higgs couplings data as shown in Fig. 8. The green, yellow, and red regions (lines) are the 68%, 95%, and 99% regions (contours), respectively. the individual panels, the constraints on the parameter space from h → τ µ and t → ch have some overlap at values of ρ = 0 and β − α = π/2. The shape of the posterior distribution in the middle panel where the constraints from the two decays are combined is determined by the BR(t → ch) measurements, since they are more accurate than the corresponding measurements of BR(h → τ µ) as is evident from Table 2. We now proceed with the combination of these flavour-violating constraints with those from the flavour conserving Higgs couplings data, discussed in Appendix B. The 2D marginalized posterior distributions are presented in Fig. 10. The 68%, 95%, and 99% contours from Fig. 8 are superimposed with dashed lines. Comparing the solid regions with the dashed lines, it is clear that the measurements of the t → ch and h → τ µ decays modify the constraints coming from the Higgs effective couplings significantly, reducing greatly the allowed degree of flavour violation encoded by the parameter ρ away from the alignment limit β − α = π/2.
The combination of constraints from Higgs data with the constraints from low-energy flavour observables that will be discussed in Appendix D is shown in Fig. 2.

D Constraints from low-energy flavour observables
The inclusive decay b → sγ is a very well measured flavour-violating process that is known to significantly constrain the parameter space of models that allow for sizable flavourviolating couplings of the Higgs boson(s) to the second and third generation quarks. Hence, it is important to understand how the parameter space of the model under consideration is affected by the measurement of this inclusive decay. The HFAG [126] average for the several measurements of the branching ratio of b → sγ is: while the SM prediction at next-to-next-to-leading-order (NNLO) stands at [127,128]: having a cutoff on the photon energy at 1.6 GeV for both the average of the measurements and the theoretical prediction. 2HDMs generically predict a contribution to b → sγ from one loop charged and neutral Higgs boson exchange. In the framework we use, one loop neutral Higgs diagrams vanish since d = 0. The charged Higgs loops will bring a new physics effect in the Wilson coefficients of the operators.
The leading-order modification to the corresponding Wilson coefficients C 7 and C 8 are given by −5y 2 j + 8y j − 3 + (6y j − 4) ln y j (y j − 1) 3 , −8y 3 j + 3y 2 j + 12y j − 7 + (18y 2 j − 12y j ) ln y j (y j − 1) 4 , The theoretical computation of the BR(b → sγ) is quite involved in the SM. In the literature, there exist extensive studies of this decay within the Type II 2HDM taking into account NNLO corrections [127]. However, such a computation is not available for a generic Type III 2HDM. It is possible to write the branching ratio in terms of modifications to the SM Wilson coefficients assuming these NP contributions are real as is the case in this study. This branching ratio with the SM contributions computed at NNLO along with the modifications to the Wilson coefficients generated by NP, δC LO 7 and δC LO 8 , calculated at the LO, can be written as: To scale the LO contribution from NP to NNLO we extract a k-factor from the branching ratios at different orders in Type II 2HDM. This factor is defined as: The k-factor does not have a significant tan β dependence in the region of parameter space we are interested in and has a weak dependence at the level of few % on the mass of the charged Higgs boson. The k-factor is fit to a polynomial function of the charged Higgs boson mass only and we assume that this functional dependence represents the k-factor for the model we consider to a good approximation, since in the limit ρ → 0 the flavour off-diagonal couplings vanish rendering the relevant couplings of the Type II 2HDM kind. The k-factor, expressed as a fourth order polynomial in m H ± , is given by with m H ± expressed in units of TeV. It should be noted here that this k-factor is only valid in the parameter space we consider in this study, specifically on the range of tan β and m H ± considered here. We also point out that while BR(b → sγ) has a sizable tan β dependence at tan β < 2 in the Type II THDM, this does not translate into a tan β dependence of the k-factor that we extract. In Fig. 11 we show the constraints on the parameter space from the branching ratio of b → sγ at the 68%, 95%, and 99% level in green, yellow and red, respectively. The

D.2 B → τ ν
In the SM, the decay B → τ ν occurs through a tree-level exchange of a W boson. In a 2HDM, this decay can also be mediated at tree-level by a charged Higgs boson, potentially giving rise to strong constraints on the parameter space. The present HFAG [126] average for the several measurements of the branching ratio is The SM prediction for the branching ratio from UTfit [129] is rather consistent with this measurement and reads BR(B → τ ν) SM = (0.807 ± 0.061) × 10 −4 . (D.10) The branching ratio of B → τ ν including NP contributions coming from charged Higgs boson exchange is given by [16,130] BR(B → τ ν) = where, τ B is the lifetime of the B + meson, f B its decay constant, and m B its mass. C ub SM , C ub R and C ub L are the Wilson coefficients of the operators where the couplings Γ H ± ij are given in Eqs. (A.3), (A.4), and (A.6). Expanding the charged Higgs couplings, it can be shown that C ub L is proportional to m u and is, therefore, vanishingly small. C ub R , on the other hand, is proportional to m b and is given by where we drop terms proportional to m µ . In the ρ → 0 limit C ub R is approximately independent of tan β. From Eq. (D.14) we conclude that the contribution to B → τ ν for ρ = 0 and tan β ∼ O(few) (as required by the constraint from b → sγ) is relatively small even for low values of m H ± . Hence, the measurement of B → τ ν does not bring a significant constraint on the parameter space of our model. Nevertheless, we will include this measurement in the global fit of the parameters discussed in Section 3.5.

D.3 R D and R D *
Beyond the B → τ ν decay, several b → c transitions have caught the attention in the last few years. In particular, there is a ∼ 4σ tension between the SM predictions for the ratios R D and R D * , and their measurement [126]. The SM predictions for these lepton universality ratios are given by Ref. [131] (see also Ref. [132]) while the several measurements of R D and R D * [133][134][135][136][137] have been combined by the HFAG collaboration [126] to give: R exp D = 0.403 ± 0.040 (stat) ± 0.024 (syst), R exp D * = 0.310 ± 0.015 (stat) ± 0.008 (syst), (D. 17) with a correlation coefficient between the two measurements of −0.23. In our 2HDM framework there can be significant contributions to both R D and R D * . A detailed analysis in a general Type III 2HDM with non-zero non-holomorphic contributions in the quark sector, but not in the lepton sector can be found in Ref. [16]. Significantly large nonholomorphic contributions and light Higgs bosons are necessary to relieve this tension between the SM predictions and the corresponding experimental values. The contribution from charged Higgs boson exchange come through the same operators O cb L and O cb R as given in Eq. (D.12) for B → τ ν, exchanging u → c. R D and R D * receive the contributions [138,139]  given in (D.12) with the exchange u → c. As in the case of the Wilson coefficient C ub R given by Eq. (D.14), C cb R is proportional to m b and is given by However, contrary to C ub L , C cb L is no longer vanishingly small and is given by where the last term derives from the u 32 contribution to the Γ H ± b L c R coupling in (A.4). Terms proportional to m µ have been dropped from both C cb R and C cb L . Despite the much larger Wilson coefficient, C cb L , if compared to C ub L , there is a necessity for a large tan β and/or quite low values for the charged Higgs boson mass, to get large enough values to explain the discrepancy between the SM values and the experimental measurements of R D and R D * [16]. Furthermore, as shown in Ref. [140], the R D and R D * anomalies prefer oppositesign Wilson coefficients, C cb L and C cb R . In the regime of large tan β, C cb R < 0, leading to the requirement for a negative value of the ρ parameter, in such a way as to obtain a positive C cb L (see Eq. (D.20)). However, as presented in Appendix D.1, the measurement of b → sγ constrains sizable values of tan β and low values of the charged Higgs boson mass and, therefore, limits the size of the NP contribution to R D and R D * to be rather small. In the combination of all constraints presented in Section 3.5, we include these measurements but do not attempt to explain the discrepancy between the SM prediction and the experimental results, which are still in their nascency.