Flavour Changing Higgs Couplings in a Class of Two Higgs Doublet Models

We analyse various flavour changing processes like $t\to hu,hc$, $h\to \tau e,\tau\mu$ as well as hadronic decays $h\to bs,bd$, in the framework of a class of two Higgs doublet models where there are flavour changing neutral scalar currents at tree level. These models have the remarkable feature of having these flavour-violating couplings entirely determined by the CKM and PMNS matrices as well as $\tan\beta$. The flavour structure of these scalar currents results from a symmetry of the Lagrangian and therefore it is natural and stable under the renormalization group. We show that in some of the models the rates of the above flavour changing processes can reach the discovery level at the LHC at 13 TeV even taking into account the stringent bounds on low energy processes, in particular $\mu\to e\gamma$.


Introduction
The second run of the LHC, at a center of mass energy √ s = 13 TeV, will provide an important probe of flavour changing couplings of the recently discovered scalar boson h [1,2]. These couplings can contribute to rare top decays like t → hq (q = u, c) and may also lead to leptonic flavour changing decays such as h → τ ± ∓ ( = µ, e), as well as hadronic flavour-changing decays like h → bs and h → bd. These decays are strongly suppressed in the Standard Model (SM) since these couplings vanish at tree level. However, Higgs Flavour Violating Neutral Couplings (HFVNC) can arise in many extensions of the SM, including in Two Higgs Doublet Models (2HDM) [3,4]. Any extension of the SM with HFVNC has to comply with the strict experimental limits on processes mediated by flavour changing neutral currents as well as with the limits on CP violating transitions leading, for example, to electric dipole moments of quarks and leptons [5].
In this paper, we investigate the allowed strength of HFVNC in the framework of a class of 2HDM, denoted BGL models, first proposed for the quark sector [6], generalised in [7] and then extended to the leptonic sector [8]. These models have the remarkable feature of having HFVNC, but with their flavour structure entirely determined by the Cabibbo-Kobayashi-Maskawa (CKM) [9,10] and the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) [11,12] matrices, denoted V and U respectively. HFVNC have been widely addressed in the literature . The distinctive feature of BGL models is the fact that many of the most dangerous HFVNC are naturally suppressed by combinations of small mixing matrix elements and/or light fermion masses. This is achieved through the introduction of a symmetry in the Lagrangian and therefore the suppression is entirely natural. Another salient feature of BGL models is that depending on the specific model of this class, HFVNC exist either in the up or in the down sector, but not in both sectors simultaneously. Analogous considerations apply to the leptonic sector. We will pay special attention to the evaluation of the maximum rate at which these processes can occur in this class of models, without violating the stringent bounds on processes like µ → eγ. In the general 2HDM, in the so-called Higgs basis [42][43][44][45], there are three neutral scalars, which we denote H 0 , R 0 and A. The couplings of H 0 to fermions are flavour diagonal in the fermion mass eigenstate basis. On the other hand, in the general 2HDM the fields R 0 and A have HFVNC with arbitrary flavour structure. The remarkable feature of BGL models is the fact that the flavour structure of HFVNC only depends on V and U . Furthermore, the neutral scalar mass eigenstates are linear combinations of H 0 and R 0 together with the pseudoscalar neutral field A. In these models the imposed symmetry restricts the scalar potential in such a way that it cannot violate CP either explicitly or spontaneously, therefore once we perform a rotation through the angle β that takes the fields from the basis chosen by the symmetry to the Higgs basis, the charged fields and the pseudoscalar field A are already physical fields. As a result, the two other neutral physical fields are related to H 0 and R 0 through a single angle rotation. In a previous work [46], we have performed a detailed analysis of the allowed mass ranges for the new scalars under the assumption that the discovered Higgs h can be identified with H 0 . We have then shown that in some of the BGL models these masses can be in the range of a few hundred GeV and thus within reach of the LHC 13 TeV run. In this paper, we work in the general case where h is a mixture of H 0 and R 0 , controlled by an angle denoted β − α. The strength of the HFVNC of h depend crucially on tan β = v 2 /v 1 , with v i the vacuum expectation values of the scalar doublets, and cos(β − α). BGL models have many features in common with several of the implementations of the minimal flavour violation hypothesis (MFV) [19,[47][48][49]. However, BGL models have the unique feature of coming from a symmetry, and as a result they have a reduced number of free parameters. This allows for definite predictions once constraints on these parameters are imposed, taking into account the present experimental bounds.
The challenge is to answer the following question: in some of the BGL models, can one have regions, in the tan β versus α − β plane, where the HFVNC of h are such that the rare processes t → hq, h → µτ can occur at a rate consistent with discovery at LHC-13 TeV? Of course, these regions have to be consistent with the stringent constraints on all Standard Model (SM) processes associated to the Higgs production and its subsequent decays in the channels ZZ, W W , γγ, bb and ττ . Furthermore the constraints derived from low energy phenomenology have to be considered: both those obtained in [46] and those new due to the presence of H 0 -R 0 mixing. Processes such as h → bs and h → bd are probably not within reach of the LHC but become important for the physics of the future Linear Collider. In this paper, we also address the question of what BGL models could lead to the observation of such decays by the future Linear Collider.
The paper is organised as follows. In the next section we briefly review the main features of BGL models in order to set the notation. In Section 3 we analyse top flavour changing decays of the type t → hq (q = u, c) in the framework of BGL models with HFVNC in the up sector. In Section 4 we perform an analysis of flavour changing Higgs decays in BGL models. In particular, we consider neutrino models with HFVNC in the charged lepton sector giving rise to h → τ decays. Up type models are also considered, giving rise to h → bs flavour violating decays. In Section 5 we investigate the discovery regions and the existing correlations for these decays in the framework of BGL models. In section 6 we present our conclusions. The paper contains two appendices, in Appendix A, we explain how the relevant Higgs experimental data has been incorporated into the analysis, and in appendix B.3, we give details relative to the low energy flavour constraint µ → eγ.

Main Features of BGL Models
Our work is done in the context of Two Higgs Doublet Models. The quark Yukawa interactions are denoted by: where Γ i , ∆ i Π i and Σ i are matrices in flavour space. The requirement that Γ i , ∆ i lead to tree level Flavour Changing Neutral Couplings (FCNC) with strength completely controlled by the Cabibbo -Kobayashi -Maskawa matrix V , was achieved by Branco, Grimus and Lavoura (BGL) [6] by means of a symmetry imposed on the quark and scalar sector of the Lagrangian of the form: where τ = 0, π, with all other quark fields transforming trivially under the symmetry. The index j can be fixed as either 1, 2 or 3. Alternatively the symmetry may be chosen as: The symmetry given by Eq. (2) leads to Higgs FCNC in the down sector, whereas the symmetry specified by Eq. (3) leads to Higgs FCNC in the up sector. These two alternative choices of symmetry combined with the three possible ways of fixing the index j give rise to six different realizations of 2HDM with the flavour structure, in the quark sector, controlled by the V matrix. We call up-type BGL models those with HFVNC in the down sector, coming from the symmetry given by Eq.
(2) and we identify each one of the three implementations by u-type, c-type or t-type depending on the value of the index j, respectively 1, 2 or 3. Likewise for the down-type models. In the leptonic sector with Dirac neutrinos there is perfect analogy with the quark sector and the corresponding symmetry applied to the leptonic fields leads to six different realizations with the strength of Higgs mediated flavour changing neutral currents controlled by the Pontecorvo -Maki -Nakagawa -Sakata matrix, U . As a result there are thirty six different implementations of BGL models. As was shown in reference [8], in the case of Majorana neutrino there are only 18 models corresponding to the neutrino types and therefore with HFVNC in the charged lepton sector. The discrete symmetry of Eqs. (2) or (3) constrains the Higgs potential to be of the form: the term in m 12 is a soft symmetry breaking term. Its introduction prevents the appearance of an would-be Goldstone boson due to an accidental continuous global symmetry of the potential, which arises when the BGL symmetry is exact. Such a potential cannot violate CP either explicitly or spontaneously. As a result the scalar and pseudoscalar neutral Higgs fields do not mix among themselves and we are left with only two important rotation angles, β and α. The angle β is such that the orthogonal matrix parametrised by this angle leads to the Higgs basis, singling out the three neutral fields: H 0 , with couplings to the quarks proportional to mass matrices, R 0 which is a scalar neutral Higgs and A which is a pseudoscalar neutral Higgs; as well as the physical charged Higgs fields H ± and the pseudo-Goldstone bosons. In BGL models the fields A and H ± are already physical Higgs fields, while H 0 and R 0 may still mix. In the limit in which H 0 does not mix with R 0 , H 0 is identified with the Higgs field h recently discovered by ATLAS [1] and CMS [2]. In this limit this field does not mediate tree level flavour changes and α is defined in such a way that the mixing angle between these fields, (β − α), acquires the value π/2. In fact expanding the neutral scalar fields around their vacuum expectation values [50] in the form φ 0 The angle β is determined by tan β ≡ v 2 /v 1 ; in the following we use the shorthand notation tan β ≡ t β , cos(β − α) ≡ c βα and sin(β − α) ≡ s βα . The general form for the Yukawa couplings of 2HDM written in terms of quark mass eigenstates and the scalar fields in the Higgs basis is given by: where γ L and γ R are the left-handed and right-handed chirality projectors, respectively, and D d and D u are the diagonal mass matrices for down and up quarks respectively. This equation defines the matrices N d and N u which give the strength and the flavour structure of FCNC and are also involved in the couplings of the charged Higgs fields. In general 2HDM, still in the Higgs basis, the flavour structure of the quark sector is fully determined by the quark masses, the V matrix and the two matrices N d and N u . It is worth emphasising the high predictive power of the general 2HDM, as can be seen from Eq. (6). Let us assume that a pair of charged Higgs H ± and the three neutral scalars (in general the physical neutral scalars are combinations of H 0 , R 0 and A) are discovered. The couplings of H ± to quarks can be readily measured from their decays. Since V in Eq. (6) stands for the CKM matrix which is known, from H ± decays one can derive N d and N u , which enables one to predict in the framework of the general 2HDM the flavour structure of the neutral scalar couplings. This would be essential to prove that the discovered neutral and charged scalar particles were part of a two Higgs doublet structure. In general 2HDM, the matrices N d and N u , are entirely arbitrary. On the contrary, BGL models have the remarkable feature of having the flavour structure of N d and N u entirely determined by fermion masses, V and the angle β with no other free parameters. This results from the symmetry introduced in the Lagrangian, in order to achieve the BGL flavour structure in a natural way. As previously emphasized, each one of the six implementations of BGL in the quark sector only has FCNC in one of the quark sectors, either up or down. In BGL up-type models the matrices N d and N u have the following simple form: where no sum in j implied. The upper index (u j ) indicates that we are considering a symmetry of the form given by Eq. (2), i.e., an up-type model with index j, thus leading to FCNC in the down-sector. Notice that all FCNC are proportional to the factor (t β + t −1 β ) multiplying products of entries involving one single row of V . The corresponding N u matrix is given by: N u is a diagonal matrix, the t β dependence is not the same for each diagonal entry. It is proportional to −t −1 β for the (jj) element and to t β for all other elements. The index j fixes the row of the V matrix which suppresses the flavour changing neutral currents. Since, for each up-type BGL model a single row of V participates in these couplings, one may choose a phase convention where all elements of N d and N u are real. For down-type models, which correspond to the symmetry given by Eq. (3), the matrices N d and N u exchange rôle, and now we have: In down-type models the flavour changing neutral currents are suppressed by the columns of the V matrix. In this paper we allow for the possibility of h being a linear combination of H 0 and R 0 which is parametrised by the angle (β − α): This mixing is constrained by data from the LHC Higgs observables (see appendix A). The quark Yukawa couplings of the Higgs field h can be denoted as: and similarly for the leptonic sector with the coefficients denoted by Y ij and Y ν ij . From Eqs. (6) and (11), we can read: More specifically, for i = j, we get the following flavour violating Yukawa couplings, for the different types of BGL models: (i) up-type u k model, with k fixed as 1 (u) or 2 (c) or 3 (t), where HFVNC arise in the down quark sector: (ii) down-type d k model, with k fixed as 1 (d) or 2 (s) or 3 (b), where HFVNC arise in the up quark sector: (iii) leptonic sector, neutrino-type, ν k model, with k fixed as 1 (ν 1 ) or 2 (ν 2 ) or 3 (ν 3 ), where HFVNC arise in the charged lepton sector: In the case of Dirac neutrinos one can write similar expressions for charged lepton type models and in this case the FCNC appear in the neutrino sector and are suppressed by the extremely small neutrino masses. In the case of models of the charged lepton type, only diagonal couplings are relevant. These couplings, as all other diagonal ones, can be extracted from equations (7) - (10), and if necessary making the corresponding changes of quarks by leptons.

Flavour changing decays of top quarks
In this section, we analyse flavour changing decays of top quarks t → hq. They can arise in down-type BGL models, where there are Higgs flavour violating neutral currents in the up sector. According to Eqs. (13) and (15), the couplings of the Higgs particle h with a top t and a light up-type quark u or c, in a model of type d ρ , can be written as One can then evaluate the corresponding t → hq decay rate. As previously mentioned, there are three types of models of this class, d ρ , depending on the column of the V matrix which suppresses the flavour changing currents. The result is Note that, apart from the global factor c 2 βα (t β + t −1 β ) 2 , every other factor in Eq. (18) is fixed once we choose the specific down-type model d ρ and the decay channel t → hc or t → hu. Therefore, for a given model, t → hq processes constrain the factor Table 1 we enumerate the decay channels and the models according  [51] is the Cabibbo angle [9] or Wolfenstein main expansion parameter [52].
It is clear that the most interesting models for t → hc are the s and b models, where the suppression is only at the λ 4 level, compared to the d model which has a strong suppression for the same decay at the λ 8 level. The d model has the curiosity that the suppression is higher for t → hc than for t → hu, unlike in the other two models. The branching ratio for t → hq in the d ρ type model is where Using the top quark pole mass m t = 173.3 GeV [51], m h = 125.0 GeV and M W = 80.385 GeV, one obtains f (x h , y W ) = 0.1306 . Considering then the upper bounds 0.79% from the ATLAS [53] and 0.56% from the CMS [54,55] collaborations, we obtain, for b and s-type models, the following constraint Notice that for this value, perturbative unitarity constraints have to be considered (see appendix B.1).

Flavour changing Higgs decays
The most interesting BGL models with HFVNC in the leptonic sector are the ν models, where there are FCNC in the charged lepton sector. As seen in the previous section for the quark sector, there are three neutrino-type BGL models, depending on the column of the U matrix which enters the FCNC in the leptonic sector. Using a notation analogous to the one of the quark sector and considering Eq. (16) for the h coupling to µ and τ , we have and the decay rate: Notice, again, the appearance of the same factor c 2 Table 2 lists the PMNS mixing matrix factors for the different νtype models. Table 2: U factors entering Eq. (23) for the different ν -type models; estimates, e.g. 1/9, 1/36, correspond to a tri-bimaximal U (except, of course, for |U e3 |); analogous information for h → eµ and h → eτ decays is provided.
The first direct search for lepton-flavour-violating decays of the observed Higgs boson performed by the CMS collaboration [56], led to the observation of a slight excess of signal events with a significance of 2.4 standard deviations. The best fit value is: which sets a constraint on the branching fraction Br(h → µτ + τμ) < 1.51% at the 95% confidence level. The ATLAS collaboration has presented a result based on hadronic   (14): Once again, it should be emphasised that once the up-type model u k is chosen, the strength of the flavour changing couplings only depends on the combination c βα (t β +t −1 β ) together with the down quark masses and V factors which are already known. The decay rate of h to pairs of quarks q i q j (i = j) is We thus have Assuming that Γ h Γ [SM] h , we can make the following estimate where Br SM (h → bb) = 0.578. The relevant CKM-related factors for h → bs and h → bd in all three u k BGL models are given in Table 3. We thus have, to a good approximation: • in models c and t, • in model u, • in all u, c and t models, We stress that, a priori, in models where there is no h → µτ constraint, one can reach values for Br(h → bs + sb) not far from 10 −1 . This can happen in charged lepton models of the charm and top types with c βα (t β + t −1 β ) ranging from 5 to 10. Again, see appendix B.1 for perturbative unitarity constraints on these values of c βα (t β + t −1 β ).

Correlations among Observables
One of the most interesting aspects of flavour violation in BGL models is the fact that, in this framework, it is possible to establish clear correlations between various flavour violating processes. As previously emphasized, one of the key features of the BGL models analysed in this paper is the presence of flavour changing neutral currents at tree level, but naturally suppressed by entries of the CKM and/or PMNS mixing matrices. Apart from these mixing matrices, in these models FCNC only depend on the values of tan β, cos(β − α) and fermion masses. However, the level of suppression depends on the specific BGL model, which in turn implies that the correlations differ from model to model. In this paper we analyse the tree level flavour violating decays involving the Higgs boson already discovered at the LHC, which were listed in the previous section. It should be pointed out that the analysis has to take into consideration the flavour conserving Higgs constraints already obtained from Run 1 of the LHC. In particular one has to comply with the measured signal strengths for the five relevant decay channels measured, to wit: h → W + W − , h → ZZ, h → bb, h → ττ , and h → γγ. They involve the flavour conserving couplings of the Higgs and put constraints in the t β vs. α − β available space. A detailed explanation of observables and the constraints is given in appendix A. An extended analysis of the phenomenology of the models under consideration, of the type presented in [46], is beyond the scope of this paper. However we also take into consideration the current bounds from the low energy processes B d −B d mixing, B s −B s mixing, K 0 −K 0 mixing and D 0 −D 0 mixing. The experimental limits on these processes can be, in principle, easily translated into limits on the combination c 2 βα (t β + t −1 β ) 2 which is relevant for our models, as will be shown in what follows. As pointed out before, there are several different types of BGL models and depending on the model under consideration we have FCNC in the down sector or in the up sector, but never in both sectors at the same time. We can thus analyse different kinds of correlations among the different HFVNC observables considered in the previous sections: • within the same quark sector: Values of the ratio in Eq. (33) for different models are shown in table 4. From Table 4 one can read that, in models of types s and b, it turns out that Br (dγ ) (t → hc) ≥ Br (dγ ) (t → hu) while in models s we have the more exotic relation Br (dγ ) (t → hc) ≤ Br (dγ ) (t → hu). The correlated allowed ranges for the branching ratios of these rare decay modes are shown in figure 1. The correlations -the lines -are associated with each particular model; for example, the purple line labeled by the letters b and τ is the correlation among the branching ratios Br (b,τ ) (t → hu) and Br (b,τ ) (t → hc), where the subscript (b, τ ) identifies completely the model -both the quark and the lepton type -. In that particular model, Br (b,τ ) (t → hc) can reach values up to a few times 10 −2 , that is up to the actual experimental upper bounds. In models (b, e) and (b, µ) -in yellow and green respectively -the correlations among Br(t → hu) and Br(t → hc) overlap along the same line with model (b, τ ), because all three models share the same V factors. The only difference between (b, τ ) and (b, µ) models in figure  1, is in the maximum allowed value of the factor c βα (t β + t −1 β ), whose origin is in the different predictions for the processes involving flavour conserving leptonic processes like h → ττ -as considered in appendix A -. Without taking into account flavour changing low energy constraints, it is clear that the models (b, τ ), (d, e), (d, µ), (s, e) and (s, µ), are the most promising models to discover new physics either in t → hc or t → hu. These models have flavour changing couplings in the up sector, therefore the Higgs coupling to uc could generate D 0 -D 0 mixing [64]. To avoid too large D 0 -D 0 mixing induced by tree level Higgs exchange, one can naively obtain the upper bounds for c βα (t β + t −1 β ) collected in table 5. Table 5: Naive D 0 -D 0 constraint on c βα (t β + t −1 β ) from tree level h exchange.
The potential effects of these constraints are represented in figure 1 with exclusion horizontal dashed lines for the corresponding models. The potential constraint in b-type models is irrelevant 5 , while this constraint could be more relevant in s and d models. Nevertheless, in these models we do not have just the 125 GeV Higgs boson, but in addition another scalar H and a pseudoscalar A. It is well known that the scalar and pseudoscalar contributions to D 0 -D 0 mixing cancel exactly in the limit of degenerate masses [65]. Note that the contributions to the oblique parameters [66] do also cancel in the limit of degenerate masses and no mixing between the standard Higgs and the other neutral scalars [67]. These considerations imply that one cannot translate into direct constraints the naive requirements on h couplings, since they can be relaxed in the presence of the other Higgses H and A. Although it is not within the scope of this paper to perform a complete analysis including the effects of the additional scalars, we illustrate how these cancellations operate in the case of meson mixing constraints, in Appendix B.2, and for the constraints coming from µ → eγ, where the cancellations are not so evident, in Appendix B.3. In the following, potential low energy flavour changing constraints are shown in the figures in the same fashion as in figure 1. It is important to keep in mind that they are indicative of which models are under pressure or else safer from that point of view.
The correlations in figure 1 correspond to the models of type down quark -charged lepton, (d i , j ), where FCNC are present in the up quark sector and in the neutrino sector. In these models, FCNC are proportional to neutrino masses and thus there are no flavour changing constraints coming from the leptonic sector. The constraints from the Higgs signals involve the diagonal couplings which do change and were taken into account as explained in appendix A. When down quark-neutrino type models are considered, |c βα (t β + t −1 β )| is also constrained by µ → eγ, τ → eγ, τ → µγ and other flavour violating processes. It can be shown that, in any ν i model, µ → eγ is the most constraining process as far as c βα (t β + t −1 β ) is concerned. We address in more detail the µ → eγ restrictions in section 5.3 and in appendix B.3.

Correlations in Up-Charged lepton models: h → bs versus h → bd
In order to have h → bq decays at a significant rate, we have to consider up-type models, u k , where FCNC occur in the down sector. In this case, following Eqs. (28)- (29), the correlations among the Higgs flavour changing decays to down quarks are given by The values of the ratio in Eq. (34) for the different up-type models are given in table  6; the correlations are represented in figure 2. The correlations in figure 2 follow from the full data analysis in appendix A, including the necessary study of the h total decay width in these BGL models. We can   table 7. Examples of models where low energy contraints can be relevant are models (u, e) and (u, µ). Once again, we must stress that the presence of other contributions in these BGL models can relax these low energy constraints. We include them here for the sake of completeness. Table 7: Note, however, that since the values in Table 7 are near 1, one cannot go too far above the dashed lines in figure 2 without taking into account perturbative unitarity (see appendix B.1).

Correlations in neutrino models
In neutrino models, we have flavour changing Higgs interactions in the leptonic sector giving rise to the interesting processes h → µ ± τ ∓ , e ± τ ∓ , e ± µ ∓ . The corresponding couplings are proportional to one of the lepton masses, therefore the transitions including a τ lepton are at least more probable by a factor (m τ /m µ ) 2 . We will concentrate in these transitions, containing a τ lepton, even if experimentally the µe channel is very interesting. These transitions are also proportional to c 2 βα (t β + t −1 β ) 2 , like all Higgs induced flavour changing transitions in these models; therefore, in down-type models, we will have perfectly defined correlations between h → µτ and t → hq; in up-type models, we will have correlations among h → µτ and h → bq. At present, as already mentioned, evidence from the CMS collaboration [56] points towards the possible observation of the decay h → µτ , which would definitely be "beyond the SM physics". These predictions could be checked by looking at the correlations with the channels t → hq for down-type BGL models, and with h → bq in up-type BGL models.
In BGL models of (d γ , ν σ ) type, the correlations t → hq versus h → µτ + τμ, following Eqs. (19)- (23), are controlled by Notice that these correlations are fixed by CKM and PMNS matrix elements. Nevertheless, there is also the ratio of the total width of the Higgs in BGL models versus the SM value. This ratio makes the correlation to depart from strict linearity depending on c βα and t β . In figure 3, we first show the plot where only the range of variation is displayed, that is the strict linear relation. In this plot, figure 3, one can see the effect of the upper bound on Br (νσ) (h → µτ + τμ). In particular, in models (b, ν 1 ), the maximum value that Br (b,ν 1 ) (t → hc) may reach is a few times 10 −3 . This value is smaller than the maximum allowed value in (b, τ ) models, presented in figure 1. As usual, we have included -with dashed lines -the naive constraints coming from the individual Higgs contribution to low energy flavour changing hadronic processes. In BGL models of (u k , ν σ ) type, we have similar expressions for the correlation among h → bq and h → µτ + τμ decays. The corresponding plot is shown in figure  4. We can observe again the effects of the measurement in the h → µτ channel. Nevertheless, as one can see, h → bs + sb branching ratios can still have values above the 10 −2 level.
Although Figures 3 and 4 show h → µτ decays, the values corresponding to h → eτ decays follow from the PMNS factors in Table 2. Notice that for h → eµ decays, an additional suppression factor (m µ /m τ ) 2 3.5 × 10 −3 is involved. It is important to stress that perturbative unitarity will not impose any further constraint on Figures 3 and 4 because Eq. (25) is at work. Several authors have noticed that µ → eγ constrains very severely the coupling h → µe via the mass unsuppressed two-loop Barr-Zee diagrams [68]. Since in BGL models all leptonic flavour changing Higgs couplings are fixed by U , masses and c βα (t β + t −1 β ), it is clear that the µ → eγ bound will translate into an important constraint on c βα (t β + t −1 β ), that has to be incorporated to the global analysis. However, in these two-loop diagrams -as in the case of the different neutral meson mixings -, not only the Higgs h can be exchanged, but also the other scalar H and pseudoscalar A will enter with the possibility to produce destructive interference between the different contributions. As detailed in appendix B.3, we have made a global analysis fully incorporating the different contributions to µ → eγ (in BGL models) in two different cases: (1) with varying free values of the masses m H and m A below 1 TeV, and (2) by imposing m H − m A ≤ 50 GeV and varying m H below 1 TeV: in both scenarios, oblique corrections can be maintained under control. To illustrate how these cancellations operate in µ → eγ, we represent the correlation among h → µτ + τμ and t → hc in models (s, ν 3 ) and (s, ν 1 ). In figure 5 we show this correlation in the full analysis, first without including the µ → eγ constraint, figure 5(a), and then, in figure  5(b), when we introduce the µ → eγ constraint as mentioned in scenario (1), that is with free values of m H and m A below 1 TeV. As figure 5 shows, the region of variation of the correlation remains essentially the same, meaning that there are cancellations at work, implying that in this kind of 2HDM, one cannot forget about the additional Higgses in order to impose the low energy constraints. Considering instead scenario (2), i.e. taking m H − m A ≤ 50 GeV and varying m H below 1 TeV, the corresponding plot is shown in figure 6.
It is then clear that if we include in the analysis the complete two-loop Bar-Zee contribution to µ → eγ, we can conclude from the observed changes and the actual level of precision, that effects are not yet relevant in the majority of BGL models. We have illustrated this result with down-type models, but the same happens in up-type models; therefore, h → bs correlations with h → µτ remain essentially unchanged without the inclusion of the µ → eγ constraint.

Conclusions
We analyse flavour changing scalar couplings in the framework of a class of two Higgs Doublet models where these couplings arise at tree level, but with their flavour structure entirely determined by the CKM and PMNS matrices. This very special structure of the scalar couplings is stable under the renormalization group, since it results from a discrete symmetry of the Lagrangian. The symmetry can be implemented in various ways, corresponding to a variety of BGL models. We pointed out that this class of models leads to New Physics with potential for being discovered at the LHC and/or at an ILC. We examine in detail rare top decays like t → hq (q = u, c) and leptonic flavour changing decays such as h → τ ± ∓ ( = µ, e), as well as hadronic flavour changing decays like h → bs and h → bd. All these decays occur in the SM only at loop level and therefore are strongly suppressed. In BGL models, the flavour violating couplings occur at tree level, but some of the most dangerous couplings are suppressed by small CKM elements. We address the question whether there are regions in some of the BGL models where these couplings are such that may lead to the discovery of rare flavour violating processes at the LHC-13TeV.
We also do a systematic study of the correlations among various observables which are an interesting distinctive feature of BGL models. In the search of these regions, we have taken into account the low energy restrictions on flavour violating processes as well as the stringent constraints on all SM processes associated to the Higgs production and subsequent decays in the various channels.
As far as the low energy flavour constraints are concerned, we agree with other authors that these cannot be imposed by assuming the dominance of the lighter Higgs contribution. This was known, in particular, in BGL models for neutral Meson-Antimeson mixing: there are important cancellations among different virtual Higgs contributions. We have illustrated this point showing how these cancellations operate in the twoloop, Higgs mediated, µ → eγ process, where these cancellations can appear in the amplitude, operating at the level of one or two orders of magnitude.
Two Higgs doublet models are one of the simplest extensions of the the SM. In general, they have a large number of free parameters and lead to scalar FCNC which have severe restrictions from low energy flavour violating processes. BGL models have the notable feature of having a small number of free parameters and achieving a natural suppression of these couplings, while at the same time allowing for the exciting scenario of having some flavour violating top and Higgs decays to occur at discovery level at the LHC-13TeV.

A Higgs signals
Besides the appearance of flavour changing couplings of the Higgs boson, as shown in equations Eqs. (12)- (13), flavour conserving couplings are modified owing to the mixing in the scalar sector, and thus a detailed analysis of the constraints on α − β and tan β that Higgs data impose is mandatory. The experimental information concerning the SM-like with mass 125 GeV discovered at the LHC is summarised in a set of signal strengths of the form where i labels the different combinations of production mechanisms and X the decay channels. Concerning production, the most relevant modes [69] are gluon-gluon fusion (ggF), vector boson fusion (VBF) and Higgstrahlung (WH & ZH); values used in the analysis are collected in where the dominant production in the 0/1 jet signal is gluon-gluon fusion -ggF:VBF ∼ 20:1 -, while in the 2 jets signal, although gluon-gluon fusion contributes the most, VBF and VH are not negligible, -ggF:VBF:VH ∼ 4:2:1 -.
where, as in µ ZZ , the 0/1 jet signal is gluon-gluon fusion dominated.
• h → ττ [74], where the 0 and 1 jet signals are dominated by gluon-gluon fusion and for the 2 jet-VBF tag, ggF and VBF production are similar.
• h →bb [75,76], Concerning the dependence of the couplings involved in the different production and decay channels, hW W and hZZ are rescaled by a factor s βα with respect to the SM for all the models. This affects VBF and VH production modes, h → W W, ZZ decays and the W -loop contribution to h → γγ, that we address later. For the couplings of h to fermions, the picture is more involved, they are modified in a model dependent manner. We remind in Tables 10, 11 and 12 the changes in htt, hbb and hτ τ with respect to the SM, where the hf f interaction is simply  Model type e, µ τ Neutrino model • the change in hτ τ only affects the branching ratio in µ τ τ i , • the change in hbb would in principle only affect the branching ratio in µ bb i ; however, production through the otherwise negligible bb → h process could be tan β or tan −1 β enhanced.
• the change in htt affects gluon-gluon production and the top loop in h → γγ decays.
The h → γγ decay deserves some attention. In the SM it is a loop induced process where virtual W and top diagrams interfere destructively. Besides the individual rescaling of both contributions, additional contributions mediated by the charged scalar H ± could also contribute. Scenarios with sizable H ± contributions to h → γγ require a specific analysis that is beyond the scope of this work. A regime with heavy H ± bosons can always be considered where this approximation is sound.
With all the ingredients in place, namely (i) the experimental constraints and (ii) the model predictions (simply expressed in terms of the different rescalings of SM couplings), a standard analysis of the {tan β, α − β} parameter space can be built. As an illustration of the effect of imposing agreement with the set of constraints on flavour diagonal Higgs couplings, we show allowed regions in the log 10 (tan β) vs. α − β plane for a few models in Figure 7. Since the overall agreement of different signal strengths with the SM is good, the region around c βα = 0 is in all cases allowed. Depending then on the particular structure of the tan β dependences in the couplings, the α − β span of the allowed regions for large or small values of tan β can be anticipated. In addition it should be noticed that, in some cases, the fluctuations departing from signal strengths equal to 1, can be in fact accommodated with α − β = π/2.

B Constraints
In this appendix we address the constraints imposed on c βα (t β + t −1 β ) by (1) perturbativity requirements for the couplings in the scalar potential, (2) tree level contributions mediated by the three neutral scalars h, H, A, to the mixing amplitude M 12 in neutral meson systems and (3) two loop Barr-Zee contributions to the µ → eγ decay branching ratio involving the flavour changing interactions of all three neutral scalars, impose on c βα (t β + t −1 β ). It is important to stress that while h → µτ , t → hc, hu or h → bs, bd, depend on c βα (t β + t −1 β ) and no other unknown parameter related to the scalar sector, the constraints analysed in sections B.1, B.2 and B.3 of this appendix, do involve new parameters like the masses m H and m A .

B.1 Perturbative unitarity
Neutral scalar masses and mixings arise from the scalar potential of the model and are related to the dimensionless quartic couplings λ i . Perturbativity requirements like λ i ≤ 4π, could have some impact on the allowed values for c βα (t β + t −1 β ). Following appendix D of [77] (here λ 5 = λ 6 = λ 7 = 0, and furthermore notice in addition that in [77] all λ i are two times our corresponding λ i in Eq. (4)), Since c βα (t β + t −1 β ) = It is then clear that, for m A ∼ v, having c βα (t β + t −1 β ) ∼ O(1) does not challenge naive perturbativity requirements like λ i ≤ 4π. For much larger values of m A , however, the situation is more involved, and only a detailed analysis including all relevant parameters can gauge the precise extent of the constraints on c βα (t β + t −1 β ) as a function of additional parameters. This is beyond the scope of the present work. Further relations, similar to (43), but involving m H and m h instead of m A , lead to the same conclusion on the perturbativity requirements for the λ i 's versus the values of c βα (t β + t −1 β ). It is important to stress, however, that the presence of other constraints overrules the potential role of imposing perturbative unitarity in the scalar potential: h → µτ alone, in neutrino type models, already requires c βα (t β + t −1 β ) 1; bounds on rare top decays t → hc, hu in models of types b and s impose c βα (t β + t −1 β ) < 5; for the remaining models constraints from meson mixings (addressed in the following), can play a more relevant role.

B.2 Mesons mixings
The contribution to the meson mixing amplitudes M 12 in BGL models is, to an excellent approximation (up to terms of order m d ms , m d m b , ms m b for up models, analogously up to mu mc , With the experimental bound [79], a simple estimate only considering the h-mediated contributions and Σ(m h ) = 2.016 (neglecting the term c 2 βα f (z h ) for α − β ∼ π/2), puts the stringent bound c βα (t β + t −1 β ) < 0.031. Nevertheless, ignoring the effect of the H and A terms is not justified. In Figure 8, we plot f (m 2 t /m 2 ), g(m 2 t /m 2 ) and Σ(m) as function of m, together with Σ(m h ) as a reference 7 . Attending to Eq.(47) and  8, sizable cancellations among the different contributions can easily be at work. A simple exercise can be helpful to illustrate such scenarios: fixing m H = 500 GeV and requiring c βα (t β + t −1 β ) > 1 (largely in excess of the bound that h alone would give), one can numerically scan the range of values of m A for which the experimental bound is satisfied. For an initial scanned range [100; 1000] GeV, values of m A satisfying the requirements cover the full range. In other words, deriving bounds from h alone is too simplistic. To remedy that and incorporate nevertheless the effect of µ → eγ, we follow two steps: (1) allowed t β vs α − β regions are first obtained from a scan of M H and M A values which fulfill the experimental bound on Br(µ → eγ) for each model, (2) the general analysis of t β vs α − β in each model is then restricted to the previous regions. For the first step, two different regimes are explored, (a) m H and m A are independent, (b) m H and m A are required to be similar, namely m A ∈ [0.8; 1.2] × m H . Although the outcome of this auxiliary analysis is a reduction of the allowed t β vs. α − β regions, prospects for the different rare decays under consideration in ν models are not significantly altered.