Higgs Quark Flavor Violation: Simplified Models and Status of General Two-Higgs-Doublet Model

We study quark flavor violating interactions mediated by the Higgs boson h. We consider observables involving a third generation quark, of both the up and the down quark sectors, like h → bs and t → ch. Using an effective field theory approach we systematically list all the possible tree-level ultraviolet completions, which comprise models with vector-like quarks and/or extra scalars. We provide upper bounds on the flavor violating transitions allowed by current limits stemming from low energy processes, such as meson mixing and b → sγ. We find that scenarios with vector-like quarks always have very suppressed flavor-violating transitions, while a general Two-Higgs-Doublet Model may have a sizeable rate. To study the latter case in detail, we perform a full numerical simulation taking into account all relevant theoretical and phenomenological constraints. Our results show that BR(t → ch) [BR(h → bs)] are still allowed at the sub-percent [percent] level, which are being [may be] explored at the LHC [future colliders]. Finally, we have found that the mild mass-splitting discrepancy with respect to the SM in the Bs meson system can be accommodated in the Two-Higgs-Doublet Model. If confirmed, it yields the prediction BR(h → bs) ≃ 10−4, if the new contribution to the mass-splitting is dominated by tree-level Higgs boson exchange.

In the Standard Model (SM), neutral flavor-changing transitions are absent at tree-level. They arise at the one loop-level with various (additional) sources of suppression like, for example, small elements of the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix or the Glashow-Iliopoulos-Maiani (GIM) mechanism. It is then clear that they constitute a privileged arena in the search for physics beyond the SM. The discovery in 2012 [1,2] of a Higgs-like scalar, h in the following, opened the possibility of exploring a new domain in neutral flavor-changing transitions and a strong experimental effort has followed, targeting processes like t → ch, uh or h →t * c (u) → W −b c (u), and potentially also h → bs, bd. We generically denote these processes as Higgs Quark Flavor Violation (HQFV). Since m t v/ √ 2, with v the electroweak symmetry breaking vacuum expectation value of the SM Higgs doublet, the Higgs-top Yukawa coupling is close to 1: if new physics is present, one may expect that such large couplings also manifest in observable transitions of the top quark to up or charm quarks mediated by h. Different studies of top flavorchanging neutral decays can be found in refs. [3][4][5][6][7][8][9][10][11][12][13][14][15]. Including those done in the so called 'flavorful models' [13,14] Current experimental bounds on the branching ratios of those processes are at the 10 −3 level (see for example refs. [16][17][18][19]): BR(t → hq) < 7.9 · 10 −3 , BR(t → ch) < 2.2 · 10 −3 , BR(t → hu) < 2.4 · 10 −3 , (1.1) at 95% CL. Limits on flavor-changing couplings of the top quark to the Z boson are also quite stringent (see refs. [16,[20][21][22][23][24]): BR(t → Zc) < 2.4 · 10 −4 , BR(t → Zu) < 1.7 · 10 −4 , (1.2) at 95% CL. Similar constraints apply to BR(t → qγ, qg) (see for example refs. [25,26]). Concerning flavor-changing couplings of h to other quarks, the LHC experiments have little direct sensitivity [27], while the ILC could in principle reach subpercent sensitivity for the branching ratios of h → bs, bd [28]; in ref. [29], it was found that BR(h → bs) can be as large as 10 −1 in Two-Higgs-Doublet Models (2HDMs). Indirect constraints can also be obtained from transitions (i.e. mixing) in the different neutral meson systems, K 0 -K 0 (ds), D 0 -D 0 (cu), B 0 d -B 0 d (bd) and B 0 s -B 0 s (bs), and from rare decays like b → sγ. The pure effective field theory approach with just the Higgs boson does not cover all the possible phenomenology relevant for HQFV, as was discussed in detail in the case of HLFV in ref. [30]. Therefore, it is crucial to also analyse in detail simplified models, whose extra particles may be subject to more stringent constraints, and to outline the models with the largest possible values of HQFV. We will concentrate on transitions involving the third and second quark generations. The paper is organised as follows. In section 2 we discuss quark flavor violation in the SM and beyond using an EFT approach. We list all possible simplified models and show how the general Two-Higgs-Doublet Model -type III -is the most promising scenario for large HQFV. In section 3, we concentrate on the relevant aspects of the latter. Flavor related constraints are addressed in section 4. A numerical analysis is then presented in section 5. Additional details are covered in the appendices.

JHEP02(2020)147 2 Quark flavor violation in the SM and beyond
In this section we discuss different aspects of quark flavor violation in the SM and beyond. We define the effective Yukawa couplings of the Higgs boson to up and down quarks as Yuk ≡ −q u y u q u h −q d y d q d h + H.c., (2.1) with summation over omitted generation indices understood: q u = (u, c, t) and q d = (d, s, b) are vectors in generation space (the quark fields are in their mass bases) and y u and y d are 3 × 3 complex Yukawa coupling matrices.

Effective field theory for Higgs quark flavor violation
In the SM the quark kinetic terms at the renormalizable level read where, under SU(2) L , Q 0 = (u 0 L , d 0 L ) are the quark doublets, and u 0 R , (d 0 R ) the up-type (down-type) quark singlets. "0" superscripts correspond to fields in a weak basis while the mass eigenstate basis is unlabelled. D denotes the covariant derivative for the different SM transformations. The SM Yukawa Lagrangian for the up and down-type quarks is where Φ = (Φ + , Φ 0 ) T is the SM Higgs doublet. Electroweak symmetry is spontaneously broken by Φ = v √ 2 ( 0 1 ), with v 246 GeV, and thus L Yuk includes mass terms The effective Higgs interactions of eq. (2.1), already written in the quark mass basis, have the simple form − mq vq qh for each quark q, with mass m q . That is, at tree-level, Higgs couplings to quarks do not violate flavor in the SM. This is an accidental symmetry of the SM, like gauge coupling universality, lepton flavor/number or baryon number, and will be violated at the loop-level or via effective operators. Indeed, at one loop, in the SM BR(h → bs) ∼ 10 −7 while BR(t → ch) ∼ 10 −15 [31] (the smallness of t → ch is due to the extra GIM suppression for virtual down quarks). Beyond eq. (2.3), the lowest dimension quark flavor-changing operators involving the Higgs field appear at dimension 6. We refer to them in the following as Yukawa operators. Denoting the scale of new physics by Λ, the effective Lagrangians for up and down quarks read respectively and They are represented in figure 1. After electroweak symmetry breaking (EWSB), the diagonalization of the complete quark mass matrices is Now the effective Yukawas of the Higgs in eq. (2.1) read the Higgs Yukawa interactions are no longer diagonal at tree-level, generating HQFV. Without loss of generality, we can use the mass basis for the up-type quarks. The quark charged current interactions read

Simplified models
In this section we discuss tree-level simplified models by "opening" the Yukawa operators for up and down quarks given in eqs. (2.5) and (2.6), respectively. This means that we give the nature and quantum numbers of the possible heavy mediators that could generate the previous operators after being integrated-out. For down quarks the operators are represented in figure 2. Tables 1 (for up quarks) and 2 (for down quarks) list all the possible simplified models. We follow the same approach used for the case of Higgs lepton flavor violation in ref. [30]. We have considered 2 extra particles at most: the new particles considered in each model are given in the second column of tables 1 and 2, where S stands for scalar and F for fermion. The (SU(2) L , Y) quantum numbers are given in the third column. In the last column, the form of the contributions to C u,d is provided in terms of the masses of the new particles (m F or m F j for fermions, m S or m S j for scalars) and of generic new physics couplings: scalar quartic couplings λ, dimensionful trilinear scalar couplings µ, µ j , and Yukawa-type couplings f q , where q = u, d, Q refers to the SM field involved (in that order of preference if two SM fields are involved) or q = VLQ for an interaction term involving two new vector-like quarks. Note that the CKM matrix V enters in the expression of C d (down quarks).

The derivative operators
Besides the Yukawa operators, there are other dimension 6 operators which generate HQFV. They involve covariant derivatives, and therefore we denote them as Derivative operators.
These are plotted in figure 3 and are listed in table 3. They are related by the equations of motion (EOM) to the Yukawa operators previously considered. This implies that, for instance, for up-type quarks their contribution to HQFV will be proportional to the quark masses. It is illustrative to consider them specifically. This is because some simple UV models directly generate them, and as we will show they are very constrained by limits from flavor-changing processes involving the Z boson. Moreover, some of the particles that generated the previous Yukawa operators also generate these ones.

Operator Topology
Particles Table 3. Tree-level topologies of the Derivative operators. The Higgs interactions in the last column correspond to the effective Yukawa couplings y q , provided in eq. (2.1), with q = u (q = d) for the first (last) four rows. Z couplings are in units of y q v/m q × e/(2c W s W ), while W ones are in units of V y q v/m q × e/(2 √ 2s W ).
We show further details regarding the generation of flavor-changing neutral currents from these operators in appendix A. The key point is that the flavor-changing neutral currents appear because of a mismatch between the quantum numbers on which the covariant derivative acts for the case of the Derivative operators and for the renormalizable kinetic terms of eq. (2.2) [30]. In table 3 we also list all the possible simplified models of the Derivative operators (third column), as well as the new Z-mediated quark flavor violating, charged-current (CC) and HQFV interactions. The Higgs interactions in the last column are given in terms of the effective Yukawa couplings y q provided in eq. (2.1). Notice that Z and W -boson interactions are independent of the quark mass involved. The chirality of the quarks involved can be understood from the operators.
We have seen using the Derivative operators that in the models with VLQ, both ZQFV and HQFV are related. We can derive the relationship among both explicitly, see also JHEP02(2020)147 refs. [30,42,44]. In models with VLQ, charged current couplings read (2.10) with V L the "enlarged" CKM matrix, and V R its right-handed counterpart (which arises when VLQ which are not SU(2) L singlets are considered). V L is now a n u ×n d matrix for n u up quarks and n d down quarks, and it is not unitary. 1 The neutral current couplings read The Z flavor-changing interactions in L Z are given by the non-unitarity of the mixing matrices. Similarly, the Yukawa couplings to h read Consider for example the ht L c R coupling; in the notation of eq. 14) The first term corresponds to the Derivative operator E

(4) u
, where a VLQ singlet (triplet) is exchanged, while the second term corresponds to E 1,2 u , where a VLQ doublet is exchanged. These contributions pick up a quark mass, as would be the case if one uses EOM to transform the operators. The last contribution corresponds to topology D of the up-quark Yukawa operator, where two types of VLQ are exchanged (a doublet plus a singlet or triplet). Therefore we can estimate the contribution as In this example, the dominant term for top HQFV is the second or the last one. For bottom HQFV clearly the last term dominates unless f VLQ v < m b . Correspondingly, the deviations from 3 × 3 unitarity of the CKM mixing matrix due to the presence of VLQ are where now f Q and m F are matrices in flavor space. JHEP02(2020)147

Estimates for models with vector-like quarks
The phenomenology of VLQ models has been scrutinised in the literature [8,[40][41][42][43][44][45][46][47][48][49][50][51][52]. For example, ref. [51] addresses in some detail constraints arising from meson mixing. For models with just VLQ, since ZQFV and deviations from 3 × 3 unitarity of CC interactions are related to HQFV, one can estimate some simple upper bounds on h → bs and t → ch. In the up sector, from eq. (2.11), the leading contribution to t → Zc (which occurs at tree-level), ignoring QCD corrections, is If the mixing with heavy VLQ is suppressed compared to the top exchange, the dominant contribution in the second term is q = t and (X u R ) tt 1, and thus we are left with Then, the leading prediction for t → ch (again, tree-level, m c → 0 and no QCD corrections) is which is two orders of magnitude smaller than the current sensitivity, eq. (1.1). For b → s transitions, although a similar reasoning would lead to straightforward bounds on the allowed values of BR(h → bs) in the context of VLQ extensions of the SM, experimental input on b − s transitions from Z → bs is much poorer than low energy constraints from B s mixing, B s → µ + µ − or b → sγ transitions. For the latter, diagrams with chirality flips in the VLQ lines dominate the processes, in an analogous way to those discussed in ref. [30] for the lepton sector, further suppressing HQFV. Typical bounds on Z bs couplings from detailed studies in the literature are below the 10 −4 level; one can thus estimate a rough upper bound on BR(h → bs) in the context of VLQ extensions

JHEP02(2020)147
Equations (2.21) and (2.22) illustrate that HQFV in extensions with just VLQ, with branching ratios forced to be below the 10 −5 level by ZQFV, are much less promising than scenarios with HQFV arising from a richer scalar sector. Moreover, the VLQ generating the Yukawa operators, always generate the Derivative operators, and therefore are subject to strong constraints. Therefore, in the following we focus on the simplest scalar scenario generating the Yukawa operators: topology A, the Two-Higgs-Doublet Model.

The general (Type III) Two-Higgs-Doublet Model
In this section we introduce the general 2HDM, also known as Type III 2HDM. Reviews addressing different 2HDMs can be found in refs. [32,[53][54][55][56]. In section 3.1 we discuss the scalar potential and in section 3.2 the Yukawa couplings. Aspects relevant for Higgs flavor-changing processes are studied in section 3.2.1.

The scalar potential
In a generic basis both Higgs scalar doublets Φ 1 and Φ 2 take VEVs denoted by v 1 and v 2 , respectively. One can rotate to the Higgs basis [57][58][59] where only one linear combination of Φ 1 and Φ 2 , denoted by H 1 , has a non-vanishing VEV, equal to 246 GeV, via the transformation where the angle β defines the mixing between the two doublets, with tan β ≡ v 2 /v 1 , and the short-hand notations s x ≡ sin x and c x ≡ cos x. We will also use t x ≡ tan x. In the Higgs basis, the doublets take the form where ϕ 1 and ϕ 2 are CP-even neutral Higgs fields, A is a CP-odd neutral Higgs field, H + is a charged Higgs field, and G + and G 0 are the would-be Goldstone bosons, which provide the longitudinal polarizations of the W + and the Z gauge bosons. The most general scalar potential is given in the Higgs basis by 2 where Λ i (i = 1, 2, . . . , 7) are the quartic couplings and M 2 ij are bare mass-squared parameters. In general, Λ 5 , Λ 6 , Λ 7 and M 12 can be complex but, by redefining H 1 and H 2 , one JHEP02(2020)147 can, for example, choose Λ 5 to be real [54]. We assume, for simplicity, a CP conserving scalar sector: all the parameters in eq. (3.3) are real.
The minimisation conditions can be used to eliminate M 2 11 and M 2 12 as independent parameters. Inserting , we obtain the squared mass of the charged scalar, 5) and the mass matrix of the CP-even neutral scalars where the mass of the CP-odd scalar is Thus, in the Higgs basis, the mass eigenstates h and H are a mixture of the CP-even states The mixing in eq. (3.8) is (3.10) It will turn out useful to obtain Λ 6 by combining eqs. (3.9) and (3.10) Eq. (3.10) and eq. (3.11) determine the sign of Λ 6 in terms of β − α. In the general 2HDM t β is not a physical parameter (see ref. [60] for a complete discussion regarding the significance of t β ). On the contrary, s β−α is a physical quantity; it needs to be sufficiently close to one (i.e., in the alignment or in the decoupling limit) so h is an adequately SM-like Higgs boson, in agreement with current observations. In this limit t 2(β−α) , and thus Λ 6 , approach zero (for Λ 6 = 0 h is exactly SM-like, with m 2 h = Λ 1 v 2 , see eq. (3.9)).

The Yukawa Lagrangian
In order to have HQFV, both scalar doublets must couple to the quarks. The most general Yukawa Lagrangian in the generic scalar basis and Y u2 are completely general 3 × 3 complex Yukawa matrices (generation indices are, again, understood and omitted). The lepton sector is assumed to be SM-like. The quark mass matrices are given by Rotating the quark fields into the mass eigenstate bases u a and d a (without "0" superscripts), M Q → D q and ξ Q →ξ Q ; without loss of generality we may work, as in section 2, in a basis where M U is diagonal with real and positive elements m u , m c , and m t . Then, the Yukawa lagrangian reads where a, b = 1, 2, 3. The correspondence with the notation in ref. [35], for a generic Yukawa couplingq JHEP02(2020)147 where P R,L = (1 ± γ 5 )/2, is provided in table 4. Due to Hermiticity of the Lagrangian, In this work we are interested in HQFV involving a third family quark. Therefore we will only consider the flavor-violating (complex) couplings inξ U, D between the third and the second families, and in addition for simplicity we set the diagonal coupling of the second generation to zero, that is, The only a priori requirement placed on the entries ofξ U andξ D is that they respect perturbativity, i.e. they are smaller than 4π.

Flavor-changing Higgs processes
In eq. (3.16), L Q includes flavor-changing couplings of h tobs,sb,tc andct controlled by the off-diagonal entries ofξ U andξ D in eq. (3.18). The h → bs decay width at tree-level, where we have neglected final state masses. The t → ch decay width at tree-level reads where we have neglected the charm mass. For the conjugate processt → hc,ξ U 32 →ξ U 23 . In the analysis of section 5, scalar decays are carried out using the inbuilt routines offered by 2HDMC [61]. The 2HDMC code does not support flavor-changing processes officially but the program is designed thoughtfully to allow for these processes. Nevertheless, some slight modifications had to be made, including promoting the Yukawa entries from real to complex. Furthermore, beyond eq. (3.19), h → bs receives QCD corrections at NLO that may increase the rate by 10-20% [29]. The 2HDMC includes QCD corrections for this process, and they are turned on in the analysis of section 5.

Constraints on quartic couplings
Since the Hamiltonian has to be bounded from below, the quartic part of the scalar potential in eq. (3.3) is required to be positive for all values of the fields and all scales. Furthermore, the considered vacuum should be the global minimum of the potential [62] (one could weaken the requirement and include a sufficiently long-lived metastable local minimum). The quartic couplings are also required to be perturbative, i.e. smaller than 4π. We also require that the scattering of the different scalars at high energies, controlled by the quartic part of the potential, respects perturbative unitarity: in particular, that the eigenvalues of the tree-level 2 → 2 scattering matrix do not yield probabilities larger than 1 (see e.g. [63][64][65], one loop corrections in a restricted 2HDM have been addressed in ref. [66]).

Higgs signal strengths
A necessary ingredient in the scalar sector is, of course, a neutral scalar with properties in agreement with the 125 GeV SM-like Higgs discovered at the LHC [1,2]. We identify it with h, and thus the first requirement is m h = (125.09 ± 0.32) GeV [72]. The width is also required to satisfy Γ h < 17 MeV following the result at 2σ presented in ref. [73]. The most relevant information for the phenomenological aspects of the 125 GeV scalar is the set of signal strengths µ XY for combined production (Y) and decay (X) channels which are factorized in "production × decay" model dependent factors The relevant production modes are gluon-gluon fusion (ggF), vector boson fusion (VBF), Higgs-strahlung (W h, Zh) and associated production with top quarks (tth); the corresponding factors are

JHEP02(2020)147
The corresponding factors for the relevant decay channels are Both κ P ggF and κ BR γγ arise from one loop amplitudes: the expressions can be found, for example, in ref. [74]. For h →τ τ , since we assume for simplicity SM-like Yukawa couplings in the lepton sector, κ BR τ τ = s 2 β−α (the experimental uncertainties in that decay channel are, in any case, large).
The experimental results (values and uncertainties) from the combined ATLAS and CMS analyses of LHC Run I data [75] are given in the following matrix: The ordering for decay channels (rows) is {γγ, ZZ, W W, τ τ, bb} and for production mechanisms (columns) {ggF, VBF, W h, Zh, tth}. For the missing entries "×" there is no measurement available in ref. [75]. In addition to eq. (3.27), we also include CMS and AT-LAS data from LHC Run II on h →bb and h →τ τ in the analysis of section 5: for h →bb, we consider CMS [76] and ATLAS [77] results for VBF production while for h →τ τ we combine ggF and VBF production following ref. [78]. Notice that the analysis of Higgs signal strengths only requires the 2HDM vs. SM modifying factors in eqs.

Flavor constraints
In order to study HQFV in the general 2HDM, flavor constraints have to be included. We discuss the most relevant constraints in the following.
In the down quark sector we focus on the process h → bs: in this case, the most stringent constraints come from the |∆B| = 2 process of B 0 s -B 0 s mixing and from the |∆B| = 1 radiative decay process B → X s γ. Since in the SM all flavor-changing processes are induced by W boson exchange, both processes occur at the one loop-level. Their GIM and loop suppressions make them highly sensitive to the presence of new physics contributions: in the general 2HDM these new contributions appear at tree-level in B 0 s -B 0 s and at one loop in B → X s γ. They are discussed in the following subsections. We do not consider other processes involving final state leptons like, e.g., B s → µ + µ − : since we focus on the quark sector, assuming SM-like tree-level couplings of scalars to leptons highly suppresses new contributions to these processes.

JHEP02(2020)147
Concerning HQFV in the up quark sector, as already mentioned, we focus on t → ch: we incorporate existing bounds at the 10 −3 level on BR(t → ch) (see eq. (1.1)). One could also consider constraints arising from D 0 -D 0 mixing. However, with the Yukawa couplings considered in eq. (3.18), the contribution to D 0 -D 0 mixing involvingξ U vanishes. 3 We do not consider constraints from t → cg, cγ processes since in this scenario they only arise at one loop while existing bounds on the corresponding branching ratios are similar to the ones for t → ch, which arise instead at tree-level.

Effective operators
We use an Effective Field Theory (EFT) approach to compute flavor constraints. An effective Hamiltonian is defined as where γ ij is the anomalous dimension matrix (ADM). The solution of this Renormalization Group Evolution equation, in vector notation, is given by where the evolution operator matrixÛ (µ, µ 0 ) is computed in terms of γ ji [80] and can be found using the publicly available Mathematica code DSixTools [81], see also ref. [82]. The CP violating mixing phase is 4 In the EFT description of B 0 s -B 0 s mixing, we adopt the usual operator basis: Since, as discussed below, W mediated contributions only affect C 1 , while the new scalar contributions affect C 2 , C 2 and C 4 , we do not factor out the usual G F and (V * ts V tb ) 2 in eq. (4.7). For the B s system the matrix elements of the operators in eq. (4.6) are Non-perturbative QCD effects [83] are encoded in the bag factors B Bs i (the vacuum insertion approximation corresponds to B Bs i → 1); they are given in table 9 in appendix B, together with the decay constant f Bs and the meson mass M Bs . The primed operators of appendix C have the same matrix elements as the unprimed ones (from parity invariance of QCD).

Standard model contribution
As anticipated, in the SM there are only contributions to the O 1 operator. The dominant contribution to C 1 (see figure 4a) is with S(x) the well-known Inami-Lim function [84]. The RGE for ∆F = 2 is given in appendix C: one can read the evolution of C 1 from the matching scale µ W ∼ m W down to µ B ∼ m B , given byη B (µ B ) = 0.862. Then (4.14) The theoretical error we choose is based upon the combination of QCD errors as laid out in table II of ref. [85], where a theoretical error of 6.2% is stated. Usingη B = 0.839 (see refs. [86,87]), ∆M Bs, SM agrees with ref. [85]. We have updated our final scan and predictions with this improved quantity, which gives ∆M Bs, SM = 1.32 × 10 −11 GeV. The SM final value in eq. (4.14) is larger than the observed one, specifically, its error translates into a 1.8σ discrepancy with the SM, as alluded to in ref. [85]. The B 0 s -B 0 s mixing phase reads β s, SM = (1.82 ± 0.11) × 10 −2 rad .

Two-Higgs-Doublet Model contributions
At tree-level, the B 0 s -B 0 s mixing process is mediated by neutral scalars h, H and A, as shown in figure 4b. The contributions to the Wilson Coefficients are [35] H, A). Beyond tree-level, there are contributions from neutral and charged scalar particles from box diagrams as shown in figures 4c and 4d; for the corresponding expressions we refer to ref. [35].

Explaining the discrepancy within the Two-Higgs-Doublet Model
Before addressing the full numerical analysis of section 5, it is interesting to study the parameter space in the 2HDM that can explain the 1.8σ deviation between the observed value of ∆M Bs and the SM prediction, see table 5. Notice that the 2HDM contribution can partially cancel the SM contribution, and therefore yield a better agreement with the lower observed value. For degenerate H and A as expected from EWPT, the tree-level contributions to the Wilson coefficients in eq. (4.16) give ∆M Bs, 2HDM tree-level We have also checked that in these regions the 2HDM is able to satisfy the observed value of the B 0 s -B 0 s mixing phase.
The 2HDM explanation of the discrepancy in terms of the tree-level contribution, also implies a prediction of BR(h → bs). For degenerate H, A, and much heavier than the light Higgs, the latter contribution to meson mixing dominates in eq. (4.17). This is true unless c β−α 0, for which in any case there is no contribution to BR(h → bs). Assuming a hierarchy in the off-diagonal Yukawas (taken real), for exampleξ D 32 ξ D 23 , so that the C 2 contribution to ∆M Bs, 2HDM dominates (and the mixed C 4 contribution can be neglected) we get from eqs. (3.19) and (4.17) where we used ∆M Bs, 2HDM = ∆M Bs, obs − ∆M Bs, SM , and Γ h 4.07 · 10 −3 GeV. The prediction is identical if the other Yukawa dominates,ξ D 32 ξ D 23 , so that C 2 dominates. On the other hand, for equal Yukawasξ D 32 =ξ D 23 , the mixed C 4 contribution cannot be neglected, and there is an extra term proportional to U 44 B Bs 4 b 4 inside the denominator of eq. (4.18), so that BR(h → bs) 6.3 × 10 −5 . As the angle β − α approaches π/2 this lower limit grows. We confirm these predictions with the scatter plots shown in figure 6, where we only have the SM plus the 2HDM tree-level contributions. We therefore conclude, that, if the observed discrepancy is confirmed, if accommodated in a 2HDM with negligible contributions at loop-level, it implies a prediction of BR(h → bs) 10 −5 -10 −4 . In our numerical scan, we indeed can accommodate somewhat lower values, when new contributions from the heavy Higgses, and/or those beyond-tree-level containing the other Yukawas, are significant. Similar studies have been done in the context of SU(5) with two Higgs doublets [88].

Radiative decays: BR(B → X s γ)
In addition to B 0 s -B 0 s mixing, we are also interested in the constraints imposed by the radiative decay B → X s γ, that is the transition b → sγ at the quark level. NNLO predictions (i.e. next-to-next-to-leading order in QCD) can be found in refs. [89,90]. In the context of the 2HDM, NNLO results can be found in ref. [91]; earlier NLO predictions [92,93] are sufficient for the scope of the present work.  The basis of operators that describes this |∆B| = 1 process includes four quark currentcurrent (O 1,2 ) and penguin (O 3−6 ) operators, together with photonic (O 7 ) and gluonic (O 8 ) dipole operators (see, e.g., ref. [94]). Effective Wilson coefficients C 7, 8[eff] are usually defined such that the perturbative contribution to BR(B → X s γ) is proportional to |C 7[eff] | 2 at leading order. Expressions for the LO and NLO contributions to the Wilson coefficients (at the matching scale µ W ∼ m W ) can be found in eqs. (16) and (17) of ref. [93]. Leading order contributions involving neutral scalars can be found in ref. [29]. 5 The perturbative 5 For comparison with the notation of ref. [93], (XY * ) φ JHEP02(2020)147 b → sγ decay rate is given by The inclusiveB → X s γ decay rate is measured with photon energies E γ > 1. 6 GeV, in which case the non-perturbative contributions relating the quark level and the meson decay rates are below the 5% level [95]. Attending to the different sources of theoretical uncertainty, in order to place constraints on the 2HDM contributions, we use the perturbative quark level decay rate in eq. (4.19) with a conservative theoretical error of 10%. The corresponding SM calculation is given in table 5 and is in very good agreement with the observed value.

Parameter scan
Given the large number of parameters of our general 2HDM (9 from the potential, 12 from the Yukawas in the 2-3 plane) we use a global fit using MultiNest [96] to scan over the allowed parameter space. We also use 2HDMC (Two-Higgs-Doublet Model Calculator) [61] to perform some phenomenological calculations. We do not include the SM one loop contribution to h → bs or t → ch, but we compute the new 2HDM contributions. We then plot our results using pippi [97]. The parameters and priors scanned over are given in table 6. We use the Higgs basis. To ensure that we carry out our scan over both quadrants in the physical angle we choose −π/2 ≤ (β − α) ≤ π/2.

JHEP02(2020)147
Parameter  Table 6. Parameters scanned over. We also indicate whether the priors are flat or log. In the Yukawa sector, i, j = 2, 3, and all other couplings are zero.
We need to provide likelihood functions L (or χ 2 = −2 ln L) to scan the parameter space of the model. To ensure that the masses of the scalars are positive, as well as to impose stability of the scalar potential, we use a hard cut-off: for a calculated value O calc and lower bound B i where χ max is large enough that the scanner effectively invalidates the point. The reverse of this may be used for an upper bound. Unitarity and perturbativity are imposed by a soft cut-off where B i is the upper bound at 68% confidence (improving the guidance provided to the scanner). For observables that have been measured we use a centered distribution with the observed value at O obs and error σ The final χ 2 -like function is built from all M bounds and N observations, For B 0 s -B 0 s mixing and B → X s γ, we sum the errors of experimental and calculated values in quadrature.

Results
To start with, we show in figure 7 the experimental contributions to the total χ 2 value that we calculate in the SM limit, that is s β−α → 1 andξ U ij =ξ D ij = 0. The largest pulls JHEP02(2020)147  here come from SM Higgs decays, as expected predominantly from h → W W , due to the fact that the experimental values of some of the production channels are slightly off from the SM, see eq. (3.27). LHC Run II data [98][99][100] gives the h → W W signal strengths by production channel (as in eq. (3.27)) as (1.10 +0.21 −0.21 , 0.62 +0.36 −0.35 , 2.3 +1.2 −1.0 , 2.9 +1.9 −1.3 , 1.5 +0.6 −0.6 ). This almost halves the h → W W channel χ 2 SM-limit contribution to ∼ 7. As such, had we included LHC Run II data in our fit we would improve our χ 2 from this degree of freedom. In any case, the SM is consistent with this data at the ∼ 2σ level.
In figure 8 we show the pull from each constraint at our best fit point for the 2HDM. This occurs at heavy scalar masses (m H = m A = m H ± ) of 2450 GeV. Relative to the χ 2 SM-limit shown in figure 7, we see that the Higgs decay channels are very similar, except for the decrease in the h → γγ channel. There is a small pull from the oblique parameters. In light of the combined fit, the pull for oblique parameters is optimised at heavier masses, JHEP02(2020)147  In the top panels of figure 9, we plot log 10 (|Λ 6 |) (left) and log 10 (c β−α ) (right) versus m H . On the top-left there is a correlation between Λ 6 and m H (as expected from eq. (3.10) for a sufficiently SM-like Higgs boson, i.e., in the alignment limit s β−α → 1). The bottom panels display correlations between the extra scalars. They each obey a linear relationship imposed by the oblique parameter constraints. The size of our masses extends up to ∼ 3200 GeV due to the priors on M 22 and the perturbativity limits used on the quartic couplings.
Exploring the constraints that caused these limits, we show in figure 11 the posterior distributions of relevant flavor physics observables (the mass splitting ∆M Bs , the CP vio-JHEP02(2020)147 lating phase β s and the radiative B-decay, B → X s γ) with respect to the h → bs decay. For ∆M Bs we observe two solution regions, as expected from figure 5. In the upper region, the predicted ∆M Bs mass splitting coincides with the SM value, which is 1.8σ away from the observed value. In the lower region, the 2HDM can accommodate the observed value, and what is more interesting, this yields a lower bound BR(h → bs), at the level of 10 −5 -10 −4 (at 1σ). This lower bound coincides well with our tree-level prediction (4.18).
In figure 12 we plot the B 0 s meson mixing mass splitting and B → X s γ versus BR(t → ch). For radiative B-decays, the combinations ξ U 23 ξ U 33 m t with tops and ξ U 23 ξ D 33 m b with bottoms in the loop, enter. On the other hand, Higgs data favours somewhat large diagonal Yukawa contributions. This in turn implies some (weak) upper bounds on ξ U 23 . The upper limit on the BR(t → ch) comes from the LHC observed upper limit, 2.2 × 10 −3 (see eq. (1.1)), hence, indirect constraints are weaker. As such, there is still almost an order of magnitude of precision before we may begin exploring the allowed 2HDM region at colliders. In this case, no lower bounds have been found from our scans, these are again just from the priors. JHEP02(2020)147   It is also interesting to investigate flavor violation in the new scalar sector, that is decays involving H, A and H ± . Figure 13 displays the modulus of the relevant off-diagonal Yukawas versus BR(H → bs). Similar plots are obtained for A → bs, tc, and H → tc. It is remarkable that these flavor-changing decays can saturate the decay widths of the heavy scalars. This may be relevant for direct searches. We also note that H + → tb has the largest lower bound.

Conclusions
In this work, we have investigated quark flavor violation involving the second and third families from an effective field theory point of view. We concentrated on the interesting processes h → bs and t → ch. After outlining the possible tree-level simplified models, which involve new scalars and/or vector-like quarks, and estimating their contributions to HQFV processes, we have focused on the most promising scenario to produce large signals: a general (or Type III) 2HDM.
We carried out a comprehensive global scan of the 2HDM imposing theoretical and experimental constraints. We focused primarily on B-physics constraints coming from B 0 s meson mixing (mass splitting and CP violating phase) and the radiative decay B → X s γ, which impose the most significant restrictions on the non-diagonal Yukawa elementsξ D

23,32
andξ U 23,32 . We have also obtained that the ∼ 2σ mass-splitting discrepancy with respect to the SM in the B s meson system can be accommodated in the 2HDM at tree-level, yielding a lower bound prediction of BR(h → bs) 10 −5 -10 −4 if loop-level and heavy Higgs contributions are not significant.
The final values obtained in our full parameter scan are BR(h → bs) < 10 −3 (10 −1 ) and BR(t → ch) < 6 × 10 −4 (10 −2 ) at 1 and 2σ (lower bounds, if present, are at the level of the one-loop SM prediction). This parameter space is already accessible and can be further examined at future colliders [101]. For example, a future 100 TeV proton-proton collider is able to constrain the t → ch channel at O(10 −5 ) [102]. Beyond the two hallmark decays, JHEP02(2020)147 possibly the easiest HQFV process to observe is H + → tb due to its large production cross section and the possible large branching fraction.
Any observed (therefore sizeable) signal of quark flavor violation involving the Higgs boson would clearly point to physics beyond the SM. As we have studied in this work, the stringent limits from low energy observables imply that it would most possibly stem from a 2HDM. We have demonstrated that the allowed parameter space in the up and the down sectors allowed by current upper limits are well within reach.

Acknowledgments
We are grateful to A. Azatov for useful discussions. This work has been supported in part by the Australian Research Council through the Centre of Excellence for Particle Physics at the Terascale CE110001004. MN acknowledges support from Fundação para a Ciência e a Tecnologia (FCT, Portugal) through the projects UID/FIS/00777/2019, CERN/FIS-PAR/0004/2017, and PTDC/FIS-PAR/29436/2017. JHG acknowledges financial support from the H2020-MSCA-RISE project "InvisiblesPlus", and he thanks the Theoretical Physics Department of Fermilab, where this project was completed, for the kind hospitality. MW is supported by the Australian Research Council Future Fellowship FT140100244.

A Derivative operators for vector-like quarks
In table 7 we list the quantum numbers (2T + 1, T 3 , Y ) of the different SM (EFT) quark objects on which the covariant derivative acts in order to derive the Z, W couplings. Details on the procedure are given in ref. [30]. In Diff. we just take the difference of the pair (T 3 , Y ) EFT − (T 3 , Y ) SM which will give us the "left-over" combination of W 3 and B fields, and therefore the Z and W interactions.

B Parameter values
The SM values used for the calculation are presented in table 8 and relevant parameters for meson mixing are given in table 9. The complex CKM matrix we use in our calculation is attained from UTFit 2016 SM Fits [103], and reads

C Evolution matrix for meson mixing
We extract the RGE matrix for the Meson Mixing basis introduced in eq. (4.6) using DSixTools [82]. This matrix represents the running of the operators from µ i = m W to JHEP02(2020)147