Role of $\lambda_{hH^+H^-}$ in Higgs boson decays $h$ to $bs$ in the 2HDM

Within the Two Higgs Doublet Model (2HDM) with $\mathcal{CP}$ conservation and a softly broken $Z_2$ symmetry, we analyze the flavor changing Higgs decays $h\to bs$ ($bs$ refers jointly to the two decay channels $b\bar s$ and $\bar b s$), where $h$ is identified with the SM-like Higgs boson discovered at the LHC. We provide a comprehensive study of the decay width $\Gamma (h \to bs)$ with particular focus on the most relevant effects from the triple Higgs coupling $\lambda_{hH^+H^-}$. Furthermore, we consider all the relevant theoretical and experimental constraints to determine which predictions for the $\mathrm{BR}\left(h\to bs\right)$ are still allowed by the current data. We find that the predictions for $\mathrm{BR}\left(h\to bs\right)$ in types II and III can be several orders of magnitude smaller compared to the SM value. In contrast, in type I and IV we find that the predicted enhancements in the decay rates with respect to the SM of up to about 70% and 50%, respectively, are still allowed. We discuss how these deviations from the SM are caused by interference effects controlled by the coupling $\lambda_{hH^+H^-}$ which can be large for very heavy $H^\pm$. To better understand the role of $\lambda_{hH^+H^-}$ in the $h \to bs$ decay we derive and analyze here the analytical results for the $hbs$ one-loop effective vertex that is generated by integrating out the heavy $H^\pm$.


Introduction
After the discovery of the Higgs boson [1,2] one of the main goals is to measure its properties with the highest possible accuracy.The current measurements, within the experimental and theoretical uncertainties, are compatible with the predictions for the Higgs boson in the Standard Model (SM).Nevertheless, the triple and quartic Higgs self-couplings still remain unmeasured, leaving plenty of room for new physics in the Higgs-boson sector.In particular, flavor-changing Higgs interactions can provide an interesting scenario to search for new physics.In the SM, flavor-changing neutral currents (FCNC) only occur at the loop level, and are thus highly suppressed.Therefore, physics beyond the SM (BSM) could predict much larger FCNC compared to the SM.One example are the decays h → bs and h → bs, which we denote together as h → bs.The first studies of the SM loop-induced hbs coupling date from the 80s [3,4], and more recent works, with the current data, report a prediction of BR (h → bs) ∼ 10 −7 within the SM [5,6].This value is very far from the expected experimental sensitivity in the short-medium term [7,8].
One of the most investigated models with Higgs-extended sectors is the Two Higgs Doublet Model (2HDM) [9][10][11].The CP-conserving 2HDM predicts five physical Higgs bosons: h and H are CP-even, A is CP-odd, and H ± are two charged Higgs bosons, with the masses denoted as m h , m H , m A and m H ± , respectively.We will assume that the lightest CP-even Higgs boson h is the SM-like boson discovered with a mass of ∼ 125 GeV.The addition of the second Higgs doublet can induce FCNC at the tree level, which are highly constrained by the experimental data.It is usual to consider a Z 2 symmetry in the 2HDM to avoid them, which leads to four different Yukawa textures [10,12], providing the so-called four 2HDM types.This symmetry can be softly broken by the parameter m 2  12 , with dimension of mass squared.Nevertheless, the 2HDM scalar sector can induce new sources of FCNC interactions at the loop level in addition to those present in the SM.At one loop, the charged Higgs boson H ± mediates these FCNC interactions, which can also involve couplings among the 2HDM scalar states.
In the case of the h → bs decay, the triple Higgs coupling λ hH + H − enters relevantly the 2HDM prediction.With our convention, the Feynman rule for the hH + H − interaction is given by −ivλ hH + H − .Therefore, the triple Higgs coupling λ hH + H − will be crucial in the 2HDM prediction, as it can additionally and significantly contribute to the flavor-changing neutral current rates.The 2HDM prediction for BR (h → bs) without tree-level FCNC was analyzed in Ref. [13].This decay was also studied in the 2HDM with an arbitrary Yukawa structure in Refs.[14][15][16].These 2HDM models1 induce FCNC at the tree level and require some fine-tuning in the Yukawa parameters to avoid other constraints from flavor observables, and thus we will not consider them here.The goal of the present paper is to update the 2HDM prediction from Ref. [13] considering the current known properties of the SM-like Higgs boson h, as well as the exclusion limits from the BSM Higgs-boson searches.Furthermore, we will emphasize chiefly the role of λ hH + H − in the prediction of BR(h → bs).From previous works, it is known that λ hH + H − can be as large as O (30) in all 2HDM types while respecting all relevant theoretical and current experimental constraints [17,18].The 2HDM prediction for h → bs can exhibit a strong constructive or destructive interference effect, depending on the sign of λ hH + H − and the 2HDM Yukawa type.
In this work, we will provide a comprehensive analysis of the effects from λ hH + H − , possibly leading to large distortions with respect to the SM prediction, while still in agreement with the current constraints to the 2HDM.To analyze in more detail the role of λ hH + H − in the h → bs decay, we will also present here a computation of the effective hbs vertex based on the full one-loop calculation and assuming a heavy H ± .This study is done by means of an expansion of the one-loop decay amplitude in inverse powers of the large mass m H ± .This will allow us to analyze in detail the most prominent effects from the very heavy charged Higgs in the flavor changing decay rates.It is found that they can yield (depending on the choice of the other 2HDM parameters) a non-vanishing constant (with m H ± ) shift with respect to the SM rates in the large m H ± limit.The effective vertex presented here, therefore, summarizes the leading (possibly non-decoupling) effects from the heavy charged Higgs boson within the 2HDM.This effective vertex is calculated for all four 2HDM types.
The structure of this article is as follows.In Sect. 2 we briefly review the 2HDM and fix our notation.In Sect. 3 we discuss the 2HDM prediction for the h → bs partial width and we analyze the role of λ hH + H − in the prediction.We first focus on the simplified scenario in the alignment limit in Sect.3.1.We also provide an effective vertex description of the 2HDM effects in this decay in the alignment limit after the performance of a large m H ± series expansion in Sect.3.2.In Sect.3.3 we turn to the study outside the alignment limit.In Sect. 4 we consider all relevant theoretical and experimental constraints for the 2HDM, using a similar analysis as in [17,18], and in Sect.4.1 we discuss the possible allowed values for BR (h → bs).In Sect.4.2 we discuss the possible experimental reach of present and future colliders in the BR (h → bs), as compared to the 2HDM predictions.Finally, in Sect. 5 we present our conclusions.

The Two Higgs Doublet Model
We assume the CP conserving 2HDM.The scalar potential of this model can be written as [9][10][11]: where Φ 1 and Φ 2 denote the two SU (2) L doublets.To avoid the occurrence of FCNCs a Z 2 symmetry [10,12] is imposed on the scalar potential of the model under which the scalar fields transform as: This Z 2 , however, is softly broken by the m 2 12 term in the potential in Eq. ( 1).The extension of the Z 2 symmetry to the Yukawa sector forbids tree-level FCNCs.This results in four variants of 2HDM, depending on the Z 2 parities of the fermion fields.Tab. 1 shows the possible couplings between the Higgs doublets and the fermions in the four different 2HDM types.

u-type d-type leptons
Table 1: Allowed fermion couplings in the four types of 2HDM.
We will study the 2HDM in the "physical basis", where the input parameters of the model are chosen to be: Here tan β := v 2 /v 1 is the ratio of the two vacuum expectations values (vevs).α is the mixing angle diagonalizing the CP-even Higgs sector.
246 GeV is the SM vev.From now on we use sometimes the short-hand notation s x = sin(x), c x = cos(x).In our analysis we will identify the lightest CP-even Higgs boson, h, with the one observed at ∼ 125 GeV.Also in some occasions we will use m, a derived parameter from m 12 and tan β, defined as, In the 2HDM, the couplings of the neutral Higgs bosons to SM particles are modified w.r.t. the SM Higgs-coupling predictions due to the mixing in the Higgs sector.Therefore, it is convenient to express the couplings of the neutral scalar mass eigenstates h i (= h, H, A) normalized to the corresponding SM couplings.We introduce the coupling coefficients ξ h i V such that the couplings to the massive vector bosons are given by: where g is the SU (2) L gauge coupling, c W is the cosine of weak mixing angle, , and m W and m Z denote the masses of the W boson and the Z boson, respectively.For the CP-even boson couplings one finds that ξ h V = s β−α and ξ H V = c β−α , whereas the coupling to the CP-odd boson is ξ A V = 0.In the Yukawa sector the discrete Z 2 symmetry leads to the following Lagrangian: where the parameters ξ h,H,A f , with f = u, d, l, are defined as: Type I Type II Type III/Flipped/Y Type IV/Lepton-specific/X The particular values of ξ f for the four 2HDM types shown in Tab. 2. Notice that the couplings ξ h,H f of the CP-even Higgs boson coupling to fermions can be understood as the value of the coupling normalized to SM prediction.
The potential in Eq. ( 1) defines the interactions in the scalar sector, also involving the SM-like Higgs boson.We define the triple Higgs couplings λ hh i h j couplings such that the Feynman rules are given by −i v n! λ hh i h j , where n is the number of identical particles in the vertex.Concretely, the coupling hH + H − , playing a dominant role in this work, is given by: In the so-called alignment limit, where c β−α → 0, all h interactions present in the SM tend to their SM values at tree level.However, even in the alignment limit some BSM interactions can be non-zero, such as hH + H − or HZA.

Γ (h → bs) in the 2HDM
In this section we will analyze the partial decay width of h → bs in the 2HDM with a softly broken Z 2 symmetry.Within this model, this observable was already computed for the first time in Ref. [13], prior to the Higgs-boson discovery.In this paper we update that prediction with the knowledge that there is a SM-like Higgs boson with a mass ∼ 125 GeV.Furthermore, we will pay special attention to the effects of the triple Higgs coupling λ hH + H − on the prediction of Γ(h → bs).
In our model set-up the only non-conserving CP effect in this process will come from the CP violating phase in the CKM matrix, that it is known to be very small.Therefore, we will compute solely the decay width of the process h → bs and multiply it by 2 to obtain the complete prediction for Γ(h → bs), to be understood as the sum of Γ(h → bs) and Γ(h → bs).The analytical and numerical results shown in this work were obtained with the help of the public codes FeynArts [19], FormCalc and LoopTools [20], where the implementation of the 2HDM was made with the public code FeynRules [21].In Appendix C we provide more details about the numerical computation of this decay width.
All one-loop diagrams2 contributing to the process h → bs in the 2HDM are depicted in Fig. 1.It should be noted that this process is loop induced because the Z 2 symmetry Figure 1: One-loop diagrams that contribute to the process h → bs.The quark q inside the loops can be any u-type quark.The diagrams with an asterisk (*) do not appear in the SM.The red dot marks the λ hH + H − triple Higgs coupling.
imposed in the 2HDM forbids flavor changing neutral currents at tree-level in the scalar sector.The total amplitude of the process can be written as: where V ij are the CKM matrix elements and A i (with i = 1...18) denotes the contribution to the amplitude of each of the 18 diagrams with the same numbering as shown in Fig. 1. p 1 and p 2 denote the momenta of the b and s quarks, respectively.The explicit expressions of the 2HDM Feynman rules involved in the computation and the amplitudes A i in the 2HDM can be found in Appendix A.
It should be noted that the particular fermions that are involved in this process are the various quarks, but not leptons, running in the loops and therefore, the prediction of the amplitude in the 2HDM types I and IV will be identical.Likewise, the prediction in types II and III will be also identical (see Tab. 2).The quark q running inside the loops in Fig. 1 can be any u-type quark, accounting for a total of 18 × 3 = 54 diagrams.All these diagrams are included in the computation of this decay.However, numerically some diagrams are more important than others.First, all contributions to the amplitude are proportional to either m s P R or m b P L , where P R,L = (1 ± γ 5 )/2 are the right/left chiral projectors.Since m s ≪ m b , the contribution proportional to P L will dominate.On the other hand, concerning the u-type quark running inside the loop, the contribution from the t quark is the leading one because it has the larger Yukawa couplings to Higgs bosons and the CKM suppression is milder with respect to other diagrams since V tb ≃ 1.As said above, the computation reported here is a complete one with the inclusion of all the 54 diagrams present in the 2HDM.However, for illustrative purposes, approximate and simpler expressions considering only the leading contributions to A i described above can be also found in Appendix A.
Finally, another important contribution in the 2HDM to the h → bs decay is the one coming from diagram 6, that depends on the triple Higgs coupling λ hH + H − .In our convention, this triple Higgs coupling is dimensionless, such that the Feynman rule of the vertex hH + H − is given by −ivλ hH + H − .This coupling is a derived parameter within our setup of the 2HDM and its expression as a function of the input parameters m h , m H ± , m 12 , tan β and c β−α can be found in Eq. (10).It should be noticed that the coupling λ hH + H − does not vanish in the alignment limit, i.e. setting c β−α = 0, and therefore neither does diagram 6.We will see later that this diagram and its interference effects with the rest of the diagrams is of crucial importance to the prediction of the decay width of h → bs.
One key fact of the process h → bs is that it is very suppressed due to the unitarity of the CKM matrix.In particular, the unitarity relation that suppresses this decay is: Therefore, any contribution in the amplitudes A i that does not depend on the mass of the u-type quark inside the loop, will factor out from the amplitude and vanish because of the unitarity of the CKM matrix.In order to find these m q -independent parts of the amplitude, is crucial to perform the Passarino-Veltman reduction to scalar one-loop integrals in the amplitude.Thus, these terms in the amplitude were removed explicitly from the numerical evaluation.
Since this is a loop induced process, the full one-loop result should be free of ultraviolet divergences.The only non-divergent diagrams are diagrams 3 to 7, both included, and diagram 12, whereas the remaining diagrams are UV divergent.The total divergence in the amplitude sums up to: where ∆ UV = 1/ϵ = 2/(4 − D), D = 4 − 2ϵ is the number of dimensions, and ϵ is a small parameter close to zero.If c β−α = 0, one recovers the same divergence as in the SM, which vanishes thanks to the unitarity of the CKM matrix, in particular due to Eq. (12).Similarly, in the generic case with c β−α ̸ = 0 the terms proportional to m 2 W also vanish because of the unitarity properties of the CKM matrix.The remaining terms proportional to m 2 q vanish because the relation (ξ d − ξ u ) (ξ d ξ u + 1) = 0 always holds in a Z 2 conserving 2HDM, see for instance Tab. 2.
The partial width of the process h → bs can be finally obtained as: where N C = 3 is a color factor, λ (x, y, z) = (x − y − z) 2 − 4y 2 z 2 and M 2 = spin |M| 2 is the spin-averaged squared amplitude of Eq. (11).It should be noted that we have neglected the possible CP violating effects in this process, as we discussed at the beginning of this section.

Γ (h → bs) in the alignment limit
In this subsection, for illustrative purposes and to simplify the analysis, we will first restrict our study to the alignment limit, i.e. c β−α = 0.Under this condition, the tree-level interactions of h with fermions, gauge bosons and Goldstone bosons recover their SM value.Therefore, diagrams containing SM particles and the SM-like Higgs boson h (i.e.diagrams without asterisk in Fig. 1) are as in the SM.In consequence, in the alignment limit, the new contributions to h → bs within the 2HDM arise only from diagrams with H ± propagating in the loops.Hence, we can write the amplitude of h → bs as: where A SM is the contribution from the SM-like diagrams in the alignment limit and A H ± is the contribution from diagrams 2, 6, 14 and 17, those that feature a H ± propagator in the loop.The remaining diagrams 4, 5, 9 and 11 vanish in the alignment limit.
In addition, for comparison, we have also computed the SM decay width under the same considerations described in the previous section and we have obtained the following value: which is in agreement with the predictions in Refs.[5,6].We will use this value throughout this work to compare consistently the 2HDM prediction with the SM.One interesting question here is, for which combination of 2HDM input parameters the SM prediction is recovered, corresponding to A H ± = 0, see Eq. (15).One could expect naively that this A H ± contribution should vanish for a very heavy H ± , since it appears suppressed by a factor m −2 H ± from each H ± propagator.In fact, that is the case for the contributions from diagrams 2, 14 and 17.However, diagram 6 could exhibit some non-decoupling effects given the fact that λ hH + H − contains a term proportional to m 2 H ± .This coupling is non-zero in the alignment limit and has the following expression (see Eq. ( 10)): For a very large value of m H ± the triangle loop function in diagram 6 scales ∝ m −2 H ± (see Eqs. ( 45), ( 78), ( 79) and ( 80)), whereas the term in λ hH + H − ∝ m 2 H ± gives rise to a nonvanishing contribution, which is constant with m H ± in this heavy mass limit3 .To suppress this contribution, one should then tune the prediction of λ hH + H − via a choice for m 12 such that λ hH + H − does not contain a term that depends quadratically on m H ± .Nonetheless, it should be noted that outside the alignment limit one can still find BSM effects even with a heavy H ± since the interaction hH + G − is also proportional to m 2 H ± .This coupling can only be zero if The numerical prediction of Γ(h → bs) as a function of m H ± , in the alignment limit, is shown in Fig. 2 for tan β = 2 (top) and tan β = 10 (bottom) for 2HDM types I/IV (left) and types II/III (right).The different solid color lines indicate several choices of the m 12 parameter.The SM prediction from Eq. ( 16) corresponds to the yellow dashed horizontal line.In this figure, one can see that the lines corresponding to m 12 = 0, 100, 200, 500 GeV tend at large m H ± to a fixed value of Γ(h → bs) ̸ = Γ SM (h → bs).As we discussed above, this is because the non-decoupling contribution from diagram 6 due to the presence of λ hH + H − .That constant value of Γ (h → bs) reached for large values of m H ± is larger than the SM value in types I/IV and smaller in types II/III.This implies that there is a constructive (destructive) interference between diagram 6 and the rest of the diagrams in types I/IV (II/III) for very large values of m H ± .Hence, those negative interference effects lead to the deep dips found in the prediction of the partial width in types II/III.Conversely, regions with m > ∼ m H ± which imply λ hH + H − < ∼ 0 show smaller Γ (h → bs) in the 2HDM than in the SM in types I/IV, while in types II/III the 2HDM rates are larger than in the SM, probing these interference patterns.These different patterns among the various types can traced back to the different signs of ξ d between 2HDM types I/IV and II/III, that enter in the Yukawa interactions of H ± (see Eq. ( 6) and Tab. 2).
On the other hand, the red solid line in each plot of Fig. 2 corresponds to a value of m 12 such that the coupling λ hH + H − in Eq. ( 17) is equal to zero 4 .Therefore in that case the only remnant 2HDM contribution in A H ± comes from diagrams 2, 14 and 17.For all 2HDM types, this line reach the SM prediction for m H ± ≳ 2 TeV, just as expected.The olive line indicates the prediction for the choice In this case, since λ hH + H − does not depend on m H ± , there is not a non-decoupling effect from diagram 6 and the prediction for Γ (h → bs) also reaches the SM value for sufficiently large m H ± .
In fact, the appearing gap at large m H ± between the predictions with a fixed value of m 12 and the red line with λ hH + H − = 0, or the olive line with λ hH + H − = m 2 h /v 2 highlights the nondecoupling contribution from diagram 6 for parameter choices that lead to λ hH + H − ∝ m 2 H ± .It also confirms that the interference effects discussed above are caused by the presence  Type II/III, tan β = 2  Type II/III, tan β = 10 Full H ± (for the choices of fixed m 12 ) and A (0) of diagram 6, and consequently, the value of λ hH + H − is crucial in the 2HDM prediction of Γ (h → bs).

One-loop effective hbs vertex from the large m H ± expansion
Given the presence of possible decoupling and non-decoupling effects in Γ (h → bs) described in the previous section, here we present a large m H ± series expansion that will capture such effects.We work again in the alignment limit, i.e, we set c β−α = 0.The analytical results will be presented in the form of a one-loop effective vertex, V eff H ± , for the hbs interaction generated from integrating out the heavy mode H ± , which is valid at large m H ± .Such an effective vertex can be useful to obtain a rough estimate of the most relevant numerical contributions from H ± to the h → bs decay rates in the 2HDM.
In order to obtain this effective vertex, we have performed a series expansion of the 2HDM contribution A eff H ± in powers of a small dimensionless variable x H ± , assuming that where m EW denotes generically any mass-dimension parameter (other than m H ± ) involved in the computation, i.e. m h , m W , m q and m 12 .Furthermore, for the results in this section, we have considered only the dominant contributions to A H ± , which are globally proportional to m 2 t m b .In summary, the validity of our assumptions in this section applies in practice to the following mass hierarchy: The case m H ± ∼ m 12 will be discussed separately below.)This hierarchy, on the other hand, (as discussed in the beginning of Sect.3) implies that the result for the effective vertex will only involve the left projector P L , since the P R contribution is always proportional to m s and can therefore be safely neglected.The analytical computation presented below was performed with the help of the public code PackageX [23].
We have computed the first two contributions in this series expansion, where H ± is the leading order contribution, i.e.O(x 0 H ± ), and A H ± is the next to leading order contribution, i.e.O(x 1 H ± ).The dots in the above equation correspond to terms of O(x n H ± ), with n = 2, 3, ... In the following of this subsection we first present the analytical results of A (0) H ± in the generic case where λ hH + H − has been replaced by the expression in the alignment limit given by Eq. ( 17).Subsequently, we compare them with the analytical results of the two cases with the particular choices of λ hH + H − = 0 and λ hH + H − = m 2 h /v 2 .The result obtained for the leading order contribution is the following: This is the non-decoupling contribution (valid under the above given assumptions) discussed in the previous section that comes solely from diagram 6.With the particular values of ξ u and ξ d from Tab. 2 in Eq. ( 20), we get the following non-decoupling effects for the four 2HDM types: Hence, A H ± has different signs in types I/IV and II/III in the large m H ± limit.The result for the next to leading order contribution is the following: with r = m 2 t /m 2 h .In the above result, the effects from H ± are suppressed by m −2 H ± and come from diagrams 2, 6 and 14.Diagram 17 is proportional to m s and therefore negligible in our approximations.A more detailed derivation of the above equations can be found in Appendix B.
In summary, and for practical purposes, the effects from the H ± in the alignment limit to the h → bs decay rates can be approximated by a one-loop effective hbs vertex V eff H ± , whose value, according to Eq. ( 11), is given by: b s h where A eff H ± , as defined in Eq. ( 19), receives contributions from leading order and next to leading order in x H ± as given in Eqs. ( 20)- (22).
The approximate prediction for Γ (h → bs) with the effects of the heavy H ± described by the effective vertex in Eq. ( 23) is also shown in Fig. 2 for the cases of fixed m 12 values.Both results with V eff H ± given by the leading order contribution, A H ± (dot-dashed lines), and by the sum of the leading and next to leading contributions, A (0) H ± (dotted lines) are shown for the different fixed values of m 12 .
One can see that the approximation given by A H ± already works quite well, specially for large values of m H ± , just as expected.For example, for m 12 = 0, this approximation is very close to the full prediction in type I/IV for values of m H ± around hundreds of GeV while for types II/III the approximation works fine from around 1 TeV.As expected, for the lines with larger values of m 12 , it is required to go to larger values of m H ± in order that A (0) H ± be a good approximation of A H ± .Also, this leading order approximation can reproduce quite well the dependence of the non-decoupling contributions for different values of tan β.In fact, this leading contribution can already describe that the gap between the SM and the 2HDM non-decoupling value increases with tan β in types II/III and decreases in types I/IV.
Turning to the prediction with the effective vertex including both leading and next to leading order contributions, A (0) H ± , we also see in Fig. 2 that it approximates better the full result than considering just the leading order, in particular in the low m H ± region where the O(m 2 EW /m 2 H ± ) corrections are more relevant.Furthermore, in types II/III, it also reproduces with a good accuracy the location of the deep interference minima present in these types, particularly if it takes place for a large value of m H ± .
Finally, we comment on the effective vertex approximation for choices of m 12 that do not lead to the dependence of λ hH + H − ∝ m 2 H ± for large values of m H ± .In Fig. 2 this is realized for the particular choice of m 2 12 = m 2 H ± sin β cos β that leads to λ hH + H − = m 2 h /v 2 , and for the choice of m 2 12 = sin β cos β m 2 H ± + m 2 h /2 that leads to λ hH + H − = 0.For these parameter choices, m 2  12 is constrained to contain a term directly proportional to m 2 H ± , i.e., the required expansion condition m 12 ≪ m H ± shown in Eq. ( 19) is not satisfied.Therefore, the consistent way to perform the series expansion in these cases is to consider the particular expression for m 12 , and the resulting value of λ hH + H − , before the series expansion.As noted above, the absence of a term proportional to m 2 H ± in λ hH + H − implies complete decoupling in the prediction of Γ(h → bs) for large values of m H ± .Thus, for both scenarios, after doing the expansion, we find a vanishing leading term: which summarizes the absence of non-decoupling effects in these two cases.However, there are still non-vanishing contributions to A H ± , which are suppressed by m 2 H ± and therefore summarize the decoupling effects in these two particular cases.If λ hH + H − = 0, the result for the next to leading order contribution is: which only receives contributions from diagrams 2 and 14.On the other hand, if H ± also receives contributions from diagram 6.In this case the result is: It should be noted that the two expressions above are different from the previous result of Eq. ( 22), as expected.In that result, there are additional decoupling contributions of O m −2 H ± coming from diagram 6 that are absent if m 12 = m H ± √ sin β cos β or λ hH + H − = 0.The numerical results of the effective vertex considering the explicit expressions of A (1) H ± in the effective vertex V eff H ± are also shown in Fig. 2. The dotted red line shows the result from the approximation in the case that λ hH + H − = 0 (Eq.( 25)), while the dotted olive line shows the result of the approximation for m = m H ± (Eq. ( 26)).Since A (0) H ± = 0 in both cases, they can reproduce the expected decoupling behavior given these parameter choices for m 12 .Additionally, the convergence to the SM prediction is also well reproduced by the effective vertex approximation considering the results for A H ± for both cases.In particular, the approximation yields a very good agreement with the complete prediction when m H ± ≳ 500 GeV.

Γ (h → bs) outside the alignment limit
Once we studied the behavior of Γ(h → bs) within the alignment limit, the next step is to consider this observable outside this limit.The fact that c β−α ̸ = 0 implies that the couplings of h to the fermions, gauge bosons and Goldstone bosons are not SM-like.Consequently, deviations with respect to the SM can arise from all diagrams.The h coupling to fermions scales with ξ h f = s β−α +ξ f c β−α (see Eq. ( 6)).Regarding the h couplings to gauge and Goldstone bosons, they are now proportional to s β−α and Yukawa-type independent, and consequently they are always smaller than or equal to the SM Higgs-boson couplings.Furthermore, another interactions involving the charged Higgs boson H ± are now also present in the prediction, like hH ± G ∓ , hH ± W ∓ ∝ c β−α , and consequently they will depend on the sign of c β−α and will be potentially small if c β−α is also small.Finally, the last new effect can enter thought λ hH + H − , that is a function of c β−α given in Eq. (10).Due to these features appearing outside the alignment limit, it is not possible to split the 2HDM result for the decay amplitude into the SM and the H ± contributions as in Eq. ( 15), and neither to write the (possibly) most relevant contribution from the 2HDM by means of just an effective vertex for very large m H ± as in Eq. ( 23).Thus, we present in this section the numerical results using the full formulas for the decay amplitude and we do not employ anymore the approximate formulas with the effective vertex discussed in the previous section.Since we investigate here the analytical dependencies of Γ(h → bs) we do not apply any experimental or theoretical constraints.This will be done in the phenomenological analysis in Sect. 4.
In Fig. 3 we show the effects of c β−α ̸ = 0 on Γ(h → bs) in the plane m H ± -tan β with m = 600 GeV for the 2HDM types I/IV (left column) and II/III (right column).The prediction for c β−α = −0.1,0.0, 0.1 is shown for both cases in the top, middle and bottom rows, respectively.The dotted pink contour lines indicates where the partial decay width has the same value as in the SM.The grey contour lines in the plots correspond to constant values of the triple Higgs coupling λ hH + H − , i.e. the parameter that controls to a large extend the size of the contribution of diagram 6.
We first focus on the plots in the alignment limit as shown in the middle row of Fig. 3.In types I/IV (middle left plot), we see that for λ hH + H − < ∼ 0 ( > ∼ 0) very important destructive (constructive) interference effects are found.These interference effects are mainly dominated by the effect of diagram 6, the one containing λ hH + H − , just as we found in the previous section.In fact, there is a narrow region where the prediction of the partial width can be below 10 −12 GeV, just in the region where λ hH + H − becomes negative.If we turn now to the prediction in types II and III (middle right plot), the situation is the opposite as compared to types I and IV.Now the strong destructive effects are found when λ hH + H − > ∼ 0, while in the region where λ hH + H − < ∼ 0 they are absent.It can also be seen that the contour lines become flat for larger values of m H ± , these are the non-decoupling effects for fixed m 12 that, as discussed in the previous section, can be approximated by Eq. ( 20).In the case of types II/III it would be necessary to go to larger values of m H ± to reach this flat behavior.Furthermore, we can also see the effects of large tan β in these plots.In types I/IV, it can be observed that the prediction of Γ(h → bs) exhibits a strong dependence on tan β for lower values of this parameter, whereas the dependence is less pronounced in the region of large tan β.This is because all diagrams with a H ± running inside the loop are proportional to ξ d = cot β, so they are very suppressed in this region of large tan β.Nonetheless, in types II/III, ξ d = − tan β and these diagrams are relevant when tan β is large.In fact, some leading contributions are proportional to ξ u ξ d , and therefore in these 2HDM types we see tan β independence in the contour lines for Γ (h → bs).
We turn now to the cases where c β−α ̸ = 0. We can understand these results as distortions from the alignment limit case.In types I/IV, we see that overall the partial width is slightly larger (smaller) when c β−α = −0.1(+0.1)with respect to the alignment limit case, specially for large values of m H ± .The reason for this behavior is that the SM-like diagrams are slightly enhanced for c β−α > 0 and diminished for c β−α > 0, and thus account for these subtle differences.Now we turn to types II/III, where we see noticeable differences between the plots with c β−α = ±0.1 and the prediction in the alignment limit, contrary to what we have seen in types I/IV.As before, the contribution from diagram 6 dominates the pattern of the prediction of Γ (h → bs) for low and moderate values of tan β < ∼ 5, where the prediction is very similar to the alignment limit case.However, in the large tan β region there is a new non-negligible 2HDM effect coming from diagram 14, which can be enhanced by ξ d = − tan β.In addition, the SM-like diagrams can be also strongly enhanced by the same parameter, like for example diagram 13.Furthermore, because of the interplay between all these contributions, the deep minimum due to the negative interference can be found around tan β ≃ 20 for lower values of m H ± and c β−α = +0.1.
It should be noted that the SM prediction, corresponding to the pink dashed contours in Fig. 3 is found in both sets of Yukawa types (I and IV vs. II and III) in different regions of the parameter space.In contrast to our discussion of the alignment limit in Sect.3.1, in these regions the sum of the SM-like diagrams do not yield the SM value: their non-SM contributions are canceled by the non-SM diagrams involving the H ± .

Allowed values of BR (h → bs) in the 2HDM
After the theoretical study of the prediction of Γ (h → bs), in this section we will analyze the possible size for the branching ratio still allowed by the current theoretical and experimental constraints of the 2HDM.The restrictions on the parameter space for the 2HDM was studied recently in Ref. [17,18].Here we follow the same procedure and only briefly summarize the considered constraints: • Electroweak precision data.New physics effects on electroweak observables from extended Higgs sectors can be parametrized in terms of the oblique parameters S, T and U [24,25].The values for the oblique parameters reported by the PDG [26] are very close to zero, meaning that the new physics contribution should remain small.Within the 2HDM, the most constraining oblique parameter is T , that is very sensitive to mass splittings between the additional Higgs bosons.Therefore, in the rest of this work, we will consider m A = m H ± , because this choice keeps the T under control at the one-loop level, even outside the alignment limit.5 • Theoretical constraints.These constraints consists in tree-level perturbative unitarity and the stability of the 2HDM potential.For the tree-level unitarity, we set an upper bound on the modulus of the eigenvalues of the two-body s-wave scattering amplitude matrix between the scalars fields [28][29][30].For potential stability, we demand the 2HDM potential to be bounded from below [31,32].Additionally, we also require that the minimum of the potential is a true minimum [33].All these constraints can be translated to inequalities involving the parameters λ i from the 2HDM potential in Eq. ( 1) and can be found in Ref. [18].Typically, these theoretical constraints impose the hardest bounds on the tripe Higgs couplings size and, in particular, on λ hH + H − .
• Higgs searches at colliders.We require that the Higgs bosons predicted by the 2HDM avoid the known experimental searches for new particles at the 95% CL.We made use of the public code HiggsBounds5.9[34][35][36][37][38][39], that determines if a particular configuration of an extended-Higgs model is excluded at the 95% CL by the existing experimental bounds.The theoretical predictions required by HiggsBounds to check the experimental bounds were computed with the public code 2HDMC [40].HiggsBounds contains experimental bounds from LEP, Tevatron and the LHC, including also all relevant results from the LHC Run 2.
• Rate measurements of the 125 GeV Higgs boson.We have identified the h state from the 2HDM with the 125 GeV SM-like Higgs boson discovered by the LHC.The properties of h should be in agreement with the experimental measurements of the mass and the signal strengths of such Higgs boson.To evaluate this agreement we made use of the public code HiggsSignals2.6 [39,[41][42][43], that performs a statistical χ 2 analysis considering all the Higgs boson signal strength and mass measurements available from Tevatron and the LHC, including most relevant results from the LHC Run 2. Again, the theoretical predictions needed by HiggsSignals to compute the χ 2 fit are obtained with 2HDMC.We will consider the allowed region to be not further away than 2σ (∆χ 2 = 6.18 for a two-dimensional scan, i.e. a plane) from the SM, where χ 2 SM = 85.76 with 107 observables.It is known that the discovered Higgs boson is in a very good agreement with the SM predictions, therefore the results from HiggsSignals will be translated into allowed values of c β−α very close the alignment limit.Stronger constraints on c β−α are expected in types II, III and IV, where the h → bb and/or the h → τ τ rates are enhanced by large values of tan β (depending on the type, see Tab. 2).
• Flavor observables.The new Higgs bosons in the 2HDM can induce relevant loop corrections to several rare b-flavored meson decays.In many BSM models the most important new contributions to these B-physics observables come from the effects from diagrams mediated by charged Higgs bosons, H ± .In our analysis, we will consider the 95% CL limits from the experimental average of BR (B → X s γ) and BR (B s → µµ), as given in [26].The theoretical computation of these BRs within the 2HDM was done with the public codes SuperISO [44,45] and 2HDMC.Stronger bounds are known to be found for all types in the low tan β region with light H ± .Furthermore, in types II and III, there is a nearly tan β-independent bound on m H ± > ∼ 500 GeV.

Results of BR (h → bs)
The results in this section will be given in terms of the ratio between the BR (h → bs) prediction in the 2HDM with respect to the SM, given by: The computation of the total width of h needed to estimate BR (h → bs) 2HDM was done with the public code 2HDMC [40].For the BR (h → bs) SM prediction we used the following value of the total width: obtained with 2HDMC via setting c β−α = 0 and λ hH + H − = 0 via a tuned value of m 12 6 .This choice allows for a consistent comparison of both quantities at the same accuracy level.
In Fig. 4 we show the prediction for R 2HDM in the plane m H ± -m in the alignment limit for m2 , for the chosen value of tan β.Since we are in the alignment limit, our prediction for the total width of h does not depend on the model type and therefore the values of BR (h → bs) are equal in types I/IV and in types II/III.
First, we comment on the shape of the allowed area by our set of constraints in the four types, whose boundaries are shown as solid black lines.For all types, there is a diagonal line, close to m H ± = m, coming from the potential stability constraints.On the other hand, the boundary close to the λ hH + H − = 30 contour corresponds to tree-level unitarity.These theoretical constraints are model independent and therefore equal for the four types.The lower bound for m H ± in all types comes from B → X s γ for types I/IV and from B s → µµ for types II/III.Finally, there are some disallowed areas in types II/III for low values of m due to BSM Higgs-boson searches.In both types, the bound below m H ± ∼ 800 GeV originates from a BSM search in the pp → H, A → τ τ channel [47], while the additional bound in type III for m H ± > ∼ 800 GeV comes from a gg → A → HZ → bbll search [48].Now, we turn to the discussion on the values for R 2HDM in Fig. 4. Overall, the contour lines for R 2HDM follow straight lines with different slopes in the m H ± -m plane.This reflects the dependence of the diagram 6 amplitude on these parameters, which is the dominant contribution to the process, as we have discussed above.Furthermore, the same interference pattern as before can be observed: in types I/IV, for λ hH + H − > ∼ 0 < ∼ 0 the interference is constructive (destructive), and the opposite happens in types II/III.Since negative values for λ hH + H − are disallowed by the theoretical constraints [17,18], inside the allowed region we find an enhancement w.r.t. the SM prediction in types I/IV (R 2HDM > 1) and a decrease w.r.t. the SM (R 2HDM < 1) in types II/III.In particular, in this plane for types I/IV, R 2HDM can reach values up to 1.4 close to the left "tip" of the allowed region, around m H ± ∼ 1 TeV and m ∼ 300 GeV (or m 12 ∼ 200 GeV).In types II/III, R 2HDM can be orders of magnitude smaller than 1, just around that same region.In principle, there should be a contour line where R 2HDM is exactly zero, but numerically that is not possible, and we find values as low as R 2HDM ∼ 10 −9 , inside and outside the allowed region.As a consequence of this destructive interference pattern, it is very difficult to find values for BR (h → bs) larger than the SM prediction within the allowed region in types II/III.
To emphasize on the relevance of the triple Higgs coupling λ hH + H − , in the prediction of BR (h → bs), in Fig. 5 we show R 2HDM in the plane c β−α -λ hH + H − for tan β = 2 and m A = m H ± = 1 TeV.Since in this figure, in contrast to our previous plots, we are using as an input parameter λ hH + H − , according to Eq. ( 10), m and m 12 should be set to the respective derived values: Furthermore, we choose m H = m + 120 GeV for types I/IV, and m H = m + 50 GeV for types II/III.These particular choices are made with the goal to obtain the maximum possible distortion for R 2HDM ̸ = 1, together with a sizable allowed region, specially in types I/IV.First, we comment on the shape of the allowed regions plotted in Fig. 5.The right curved  boundary in types I/IV is present due to tree-level unitarity constraints.The curved right border in type I covering values of λ hH + H − from 0 to 20 is due to the potential stability constraints.The additional upper irregular bounds for λ hH + H − > 20 in type I originate from an analysis combining several pp → H → W W, ZZ searches [49].The "left" and "right" limits on c β−α ∼ ±0.05 in types II/III and also the "left" limit in type IV are given by the LHC rate measurements of the 125 GeV Higgs boson.The upper bound for λ hH + H − ∼ 30 in type II come from a search in the pp → H → τ τ channel [47].In the case of type III, this bound is not relevant in that region and the upper limit on λ hH + H − appears due to the tree-level unitarity conditions.In addition, the small cracks in the allowed region present in type III around λ hH + H − ∼ 30 come from searches in the channel pp → H → ZZ → 4l, 2l2q, 2l2ν with l = e, µ [50].In all types, the bound for low values of λ hH + H − are present due to the potential stability conditions.Now we turn to the predictions for R 2HDM in types I/IV, shown in the left panels of Fig. 5.For these types, the contour lines for this variable follow straight lines with a constant slope, where R 2HDM is larger for lower values of c β−α and larger values of λ hH + H − .This is in agreement with our discussion in Sect.3.3, where negative values of c β−α slightly enhance the prediction of BR (h → bs) because of the effects from the SM-like diagrams.In addition, the total decay width of h has also a small effect on R 2HDM in that direction because it decreases with c β−α in types I/IV.The predictions of R 2HDM look very similar for both types, but there is a small, hardly visible difference in the total width of h, and therefore in R 2HDM , due to the differences in the leptonic decays.Within the allowed region, the maximum values of R 2HDM are reached in both types for large (but allowed) values of λ hH + H − ∼ 20 − 30 and the lowest allowed value of c β−α .Concretely, R 2HDM can reach values up to 1.7 in type I and 1.5 in type IV.The choice for m H in types I/IV is such that the region allowed by the theoretical constraints is "shifted" to negative values of c β−α , where BR (h → bs) is larger.Larger values for R 2HDM are difficult to obtain while satisfying all considered constraints.For instance, the prediction for R 2HDM could be larger for smaller values of tan β, but that would imply heavier H ± values to avoid flavor constraints.Such large values of m H ± , together with large values of λ hH + H − , easily violate the constraints from tree-level unitarity, especially outside the alignment limit.
Finally, we comment on R 2HDM in types II/III, shown in the right column of Fig. 5.Both Yukawa types yield very similar predictions.Again, there is a tiny difference in R 2HDM between both types due to different h decay rates to leptons.In this case, the predictions for R 2HDM decrease for large values of λ hH + H − , reaching a minimum in the upper region of the plot.Again, the reason is the destructive interference governed by diagram 6 in the prediction of Γ (h → bs).The contour lines are nearly c β−α -independent, but the tilt is slightly stronger for larger values of λ hH + H − .These extremely low values of the BR (h → bs) prediction can also be found in other regions of the parameter space of types II and III, in contrast to the maximum values found in types I and IV, which are only found in the indicated parameter space.

Experimental prospects for h → bs
In this section we discuss on the future prospects to detect experimentally the decay channel h → bs.It is very challenging for the LHC, as well as the future HL-LHC, to measure any flavor-changing decays in the quark sector for the SM-like Higgs boson.These searches suffer from extremely large hadronic backgrounds and therefore, in practice, they are inaccessible for these machines, see for instance Ref. [7].However, possible future e + e − colliders provide a much cleaner environment and a priori could be able to detect this process.In particular, there is an analysis of the experimental prospects at the International Linear Collider (ILC) [8].It was found that the ILC could set an upper bound on BR (h → bs) of order 10 −3 at the 95% CL.To our knowledge, there are no further analyses of this type for other lepton colliders, such as CLIC, FCC-ee or CEPC.The possible future upper bound of O(10 −3 ) is five orders of magnitude above the SM prediction.Furthermore, in the previous section, we found that only enhancements up to 70% or 50% w.r.t. the SM (i.e.R 2HDM = 1.7 and 1.5) could be realized within 2HDM types I and IV, respectively.In the case of 2HDM types II and III, the predictions are in general below the SM prediction in almost the entire parameter space.In consequence, it appears that there is no collider projected at the medium-large term with sufficient experimental reach to measure such values of BR (h → bs), neither within the SM, nor in the 2HDM.Conversely, any signal detected in the h → bs channel at the projected future experiments would clearly point to models beyond the SM and the 2HDM with a softly broken Z 2 symmetry.
One could ask the question if the prediction for R 2HDM could be larger under some circumstances.The main limiting factor to have a larger decay width for h → bs are the flavor constraints.In particular, for types I and IV, B → X s γ sets the strongest bounds at low tan β in our analysis.If in the future this bound relaxes (due to a shift in the central experimental value or any other reason), the 2HDM could allow larger rates for h → bs.For instance, from Fig. 3 it can be seen that Γ (h → bs) ∼ 10 −9 GeV for tan β = 1 and m H ± = 1 TeV, leading to R 2HDM ∼ 2.5.Nevertheless, it does not appear likely that such bounds would relax below tan β = 1.On the other hand, large enhancements for very light H ± could also be realized coming from diagrams without λ hH + H − .However, these low values of m H ± are strongly disfavored by experimental searches, see, e.g.Ref. [51] (LEP) and Refs.[52][53][54] (LHC).Consequently, only a slightly larger prediction for R 2HDM ∼ 3 could be possible if the flavor bounds on the low tan β region become less constraining in the future.However, even such enhancements are still very far away from the expected experimental reach for this process.Therefore, the conclusion holds that if any evidence of the decay h → bs is detected in the near future, this would imply the presence of new physics beyond the SM and the CP conserving 2HDM with a softly broken Z 2 symmetry.
We would like to comment briefly on the case that the discovered 125 GeV Higgs boson is identified with the heavy CP-even boson H, implying that m h < 125 GeV and that alignment limit corresponds to c β−α → 1.This is known in the literature as "inverted" hierarchy (m h < m H ∼ 125 GeV), compared to the usual "normal" hierarchy (m h ∼ 125 GeV < m H ). We do not expect any significant impact on our conclusions for BR(H → bs) in the "inverted" scenario compared to BR(h → bs) in the "normal" one.In the alignment limit, the prediction of the decay rates for H → bs in the "inverse" hierarchy is identical to h → bs in the "normal" hierarchy but with m H → m h .Outside the alignment limit, there are some contributions in H → bs that change sign with respect to h → bs (like the vertices hH ± G ∓ vs. HH ± G ∓ ), but all of these would be small since they are proportional to s β−α ∼ 0. The similarity between h → bs and H → bs in the "normal" and "inverted" hierarchies respectively is also observed in Ref. [13].Furthermore, the fact that both CP-even Higgs bosons are light in the "inverted" scenario would imply an upper bound on the mass of m H ± below 1 TeV (see, for instance, Ref. [55]).Therefore, we expect very similar conclusions in both scenarios.In addition, one could also think that the "inverted" scenario could favor lower values for m H ± , but as discussed above, such lower values of m H ± are currently very discouraged by direct LEP/LHC searches.
Finally, we would like to highlight the differences between our work and the previous results of Ref. [13] from 2004.As discussed before, we update and complement the results of Ref. [13], which investigated the BR(h → bs) only in the 2HDM types I and II.In the present work, we have analyzed the allowed values for BR(h → bs) in the four 2HDM types considering all the latest relevant constraints, including the measured properties of the later discovered Higgs boson with a mass around 125 GeV.Additionally, we also highlight the relevant role of the triple Higgs coupling λ hH + H − in the final prediction of BR(h → bs).We discuss how this observable can exhibit non-decoupling effects induced by λ hH + H − , and we have also computed an effective vertex hbs valid for large values of m H ± in the alignment limit.

Conclusions
Within the 2HDM we have analyzed the phenomenology of the flavor-changing decays of the SM-like Higgs boson to a strange-bottom quark pair, h → bs and h → bs, which we denote together as h → bs.We have assumed the lightest CP-even Higgs boson, denoted as h, to be the SM-like Higgs discovered by the LHC.We consider a Z 2 symmetry that forbids FCNCs at the tree level, which is softly broken by the parameter m 12 .This is in contrast to other recent works that study this process within a 2HDM with tree-level flavor-changing Yukawa interactions [14][15][16], which can predict more sizable rates, but has to be fine-tuned to avoid large FCNCs disallowed by the current experiments.Hence, this work provides an updated analysis of the previous result for one-loop generated FC Higgs decays within the 2HDM from Ref. [13], published before the discovery of the Higgs boson at ∼ 125 GeV.Furthermore, a significant focus of this paper is the analysis of the 2HDM contributions mediated by the charged Higgs bosons H ± .More concretely, we emphasize the role of the triple Higgs coupling λ hH + H − , which was found to be crucial in this decay process.
We started our analysis of the h → bs decay width in the alignment limit, which implies that the 2HDM couplings of the h that are present in the SM recover their SM values.Under this assumption, we found that the amplitude for the decay process h → bs can be split into two parts.The first part corresponds to the set of diagrams as in the SM, while the second one consists of a pure 2HDM set of diagrams involving the charged Higgs boson, H ± .Among these 2HDM diagrams, the one with the triple Higgs coupling λ hH + H − stands out because its contribution can be of a similar order or even larger than the SM-like contribution.Consequently, the interference between the diagram with λ hH + H − and all other contributions turns out to be crucial in the analysis of the prediction of the decay width.Concretely, this interference is positive (negative) for λ hH + H − > 0 in the 2HDM types I and IV (types II and III).Furthermore, this diagram can possibly lead to non-decoupling effects that do not vanish even for heavy H ± due to the dependence on m H ± of the triple Higgs coupling λ hH + H − .We captured these one-loop effects of a heavy charged Higgs boson in the prediction of Γ(h → bs) by a series expansion in inverse powers of m H ± up to (m EW /m H ± ) 2 , assuming this mass to be much larger than all other occurring masses.Using this series expansion, we presented a one-loop effective vertex description of the heavy charged Higgs effects in h → bs in the alignment limit.This effective vertex approximation yields reasonably accurate predictions of the decay width for large values of m H ± above 500 GeV, and are particularly good for low values of tan β < 10.
Subsequently, we analyzed the effects outside the alignment limit, obtaining similar conclusions since the effect of moderate values of c β−α ̸ = 0 slightly modifies the h couplings to other particles.Therefore, the interference pattern between the diagram with λ hH + H − and the remaining contributions holds similarly as in the alignment limit.
We studied the possible range of BR(h → bs) allowed by all relevant theoretical and current experimental constraints on the 2HDM parameter space, based on previous analysis from Refs.[17,18].Mainly positive values of λ hH + H − are allowed by the 2HDM potential stability constraint.Consequently, due to the interference patterns found in the decay width prediction for the 2HDM Yukawa types, we find enhancements of BR (h → bs) w.r.t. the SM prediction in types I and IV and a decrease in the allowed values for the BR in types II and III.More concretely, in types I and IV, one can find enhancements of the BR prediction of around 70% and 50% w.r.t. the SM, respectively.These predictions are realized for negative values of c β−α and large values of λ hH + H − .The principal constraints that yield this maximum in the BR predictions are light Higgs-boson rate measurements at the LHC, which do not allow for a large deviation from the alignment limit, as well as the theoretical constraints on the potential, due to the large values of λ hH + H − .Regarding types II and III, the negative interference effects in these types lead to significantly suppressed BR (h → bs) by up to several orders of magnitude compared to the SM.
Finally, we discussed the prospects to observe and probe the h → bs decay experimentally.We conclude that the allowed ranges for the 2HDM prediction of BR (h → bs) are several orders of magnitude beyond the experimental reach of present and even future collider experiments.Therefore, in the hypothetical case that any signal from this h decay channel is detected experimentally, it would require from physics beyond the SM and the 2HDM with a softly broken Z 2 symmetry.
with λ hH + H − defined in Eq. ( 10) and all momenta going inwards.
The final contribution from diagrams 2, 6, 14 and 17 in the alignment limit up to order one in x H ± is: where we have neglected the m t -independent terms, because they vanish due to the unitarity of the CKM matrix.Under the assumption that m s = 0 we get: Finally, if m b is neglected with respect to m t , m h and m 12 we get:

C Numerical input parameters and considerations
The correct implementation of the unitarity of the CKM matrix is key to compute the partial width of the decay h → bs.Therefore, we parametrized the CKM matrix V CKM in the "standard parametrization" consisting in three rotation angles and a CP phase δ, yielding a CKM matrix that is unitary by construction.Furthermore, we considered the CP phase to be zero, since the CP effects in the observable h → bs are negligible.By doing this, the CKM matrix becomes real and the unitarity conditions are easier to fulfill numerically.In summary, the CKM matrix that we used in our computation is: where s 12 = 0.22500, s 23 = 0.04182, s 13 = 0.00369, according to PDG [26].In order to ensure with high numerical accuracy the unitarity of the CKM matrix, we have computed the Passarino-Veltman functions present in the amplitude (shown in the previous appendix) with the help of LoopTools [20] with quadruple precision.
In addition, we also considered the masses of the final-state quarks, m b and m s , in the MS scheme at an energy scale of the SM-like Higgs boson mass, in order to include leading contributions from QCD corrections.We used the values as given in Ref. [56]: m b (m h ) = 2.768 GeV and m s (m h ) = 0.052 GeV.The rest of the masses involved in the computation were taken from PDG [26]: m u = 2.16 × 10 −3 GeV, m c = 1.27GeV, m t = 172.5 GeV, m W = 80.377 GeV and m h = 125.25 GeV.In our calculation we employed the G F -scheme, where the EW parameters are derived from the Fermi constant G F = 1.1663788×10 −5 GeV −2 .In particular, v = ( √ 2G F ) − 1 2 ≃ 246.22 GeV and g = 2m W /v ≃ 0.65.

Figure 2 :
Figure 2: Prediction of Γ(h → bs) in the 2HDM types I and IV (left) and types II and III (right) as a function of m H ± in the alignment limit for tan β = 2 (top) and tan β = 10 (bottom) and for different values of m 12 .Dot-dashed and dotted lines correspond to the 2HDM prediction of Γ(h → bs) with the charged Higgs boson effects approximated by a large m H ± expansion given by an effective vertex V eff H ± , Eq. (23) with A H ± respectively.Dashed yellow horizontal lines indicate the SM prediction.

Figure 3 :
Figure 3: Prediction of Γ(h → bs) in the plane m H ± -tan β with m = 600 GeV for 2HDM types I/IV (left column) and types II/III (right column).c β−α is set to -0.1 (top), 0 (middle) and 0.1 (bottom).Dotted pink line corresponds to the SM prediction.Solid grey lines show contour lines for different values of λ hH + H − .

Figure 4 :
Figure 4: Prediction of R 2HDM = BR (h → bs) 2HDM /BR (h → bs) SM in the plane m-m H ± with m h = 125 GeV, m H = m + 50 GeV, m A = m H ± , tan β = 2 and c β−α = 0 for the four 2HDM types.The black contour lines indicates the boundary of the allowed region by all considered constraints.The dashed pink contour correspond to R 2HDM = 1.Solid grey lines shows contour lines with different values of λ hH + H − .

tan β = 2 ,
m h = 125 GeV, m H = m − 50 GeV and m A = m H ± .The R 2HDM prediction for types I/IV are shown in the left column of the figure, while the 2HDM types II/III predictions are shown in the right column.In the upper axis we show the corresponding value of m 12 , which is related to m by the equation m 2 12 = m2 sin β cos β = 2 5

Figure 5 :
Figure 5: Prediction of R 2HDM in the plane λ hH + H − -c β−α with m h = 125 GeV, m H ± = 1000 GeV, tan β = 2 and m H = m+120 GeV for 2HDM types I and IV and m H = m+50 GeV for 2HDM types II and III.The black contour lines indicates the boundary of the allowed region by all considered constraints.The dashed pink contour correspond to R 2HDM = 1.

A 4 = c β−α m 2 h − m 2 H ± m s P R m 2 b 2 −
m b P L ξ d m 2 s C

Table 2 :
Values for ξ f in the four Z 2 conserving 2HDM types.