Two Higgs Doublet Models and $b\to s$ exclusive decays

We computed the leading order Wilson coefficients relevant to all the exclusive $b\to s\ell^+\ell^-$ decays in the framework of Two Higgs Doublet Model (2HDM) with a softly broken $\mathbb{Z}_2$ symmetry by including the $\mathcal{O}(m_b)$ corrections. We elucidate the issue of appropriate matching between the full and the effective theory when dealing with the (pseudo-)scalar operators for which keeping the external momenta different from zero is necessary. We then make a phenomenological analysis by using the measured $\mathcal{B}(B_s\to \mu^+\mu^-)$ and $\mathcal{B}(B\to K \mu^+\mu^-)_{\mathrm{high}-q^2}$, for which the hadronic uncertainties are well controlled, and discussed their impact on various types of 2HDM. A brief discussion of the decays with $\tau$-leptons in the final state is also made.


Introduction
Physical processes driven by the flavor changing neutral currents (FCNC) are forbidden in the Standard Model (SM) at tree level. Since they occur through loops, their measurements offer a low energy window to the particle content in the loops. In other words, they do not only represent a fine test of validity of the SM, but they also offer an opportunity to look for the effects of physics (particles) beyond the SM (BSM) at low energies. The main obstacle to the accurate comparison between the SM theory and the experimental data lies in the fact that the non-perturbative QCD effects are not under full theoretical control. While the solution to non-perturbative QCD is lacking, in some situations the hadronization effects can be solved by means of numerical simulations of QCD on the lattice (LQCD). Over the past couple of decades we witnessed a huge progress in reducing the uncertainties in the LQCD results. Nowadays, an excellent theoretical control of the neutral meson mixing processes promoted those FCNC processes to viability tests of the New Physics (NP) model candidates. Besides the oscillation frequencies of the neutral meson systems, the processes based on b → s transitions received a great deal of attention in the particle physics community. While the inclusive and exclusive processes based on the penguin-induced b → sγ decay have been, and still are, a very significant constraint when building a NP model, the processes based on b → s + − received a huge attention because they allow to access another types of penguin and box diagrams. With the advent of the Large Hadron Collider (LHC) the measurement of B(B s → µ + µ − ) became possible [1] and the result appeared to be somewhat lower than predicted [2]. The spectrum of dB(B → Kµ + µ − )/dq 2 has been measured [3] too and in the range of large q 2 's it appears to be larger than predicted [4,5]. A full angular analyses of B(B → K * µ + µ − ) [3,6] and B(B s → φµ + µ − ) [7] revealed discrepancies in several observables with respect to their SM predictions [8]. Moreover, the measurement of R K = B (B → Kµ + µ − )/B (B → Ke + e − ) [9] was shown to be significantly lower than predicted (by 2.4σ) [10]. 1 Those new experimental data helped discarding several NP models and are currently used as constraints in building a NP model.
Simultaneously with the research of FCNC processes, the LHC experiments allowed observing the missing ingredient of the SM, the Higgs boson, the mass of which has been found to be m h = 125.09 (24) GeV [11]. While this was a milestone of the LHC, the pending question of hierarchy of scales remains open and a quest for physics BSM continues. One of the minimalistic approaches to building a model of physics BSM is to extend the Higgs sector by introducing an extra Higgs doublet. Phenomenology in the scenarios with two Higgs doublets appears to be very rich and the associated models are generically called the Two Higgs Doublet Models (2HDM), cf. e.g. [12][13][14]. Nowadays the experimental search of the additional Higgs bosons is one of the main goals at LHC, in particular that of the charged Higgs boson [15]. Like in the SM, introducing fermions to the 2HDM context results in a plethora of new parameters. To restrain the number of those parameters and to prevent from appearance the FCNC at tree level it is common to assume a peculiar pattern of Yukawa couplings. To test those assumptions one needs to compare many measured observables with theoretical expressions derived in SM with the extended Higgs sector. In this paper we elaborate a few lessons one can learn from the measured b → sµ + µ − processes about 2HDM with a softly broken Z 2 symmetry. In doing so we will use two observables, namely B(B s → µ + µ − ) and B(B → Kµ + µ − ) high−q 2 , which are very well measured experimentally and for which the theoretical control of the corresponding hadronic uncertainties is established by the LQCD computations [16]. For other observables the theoretical uncertainties are not as well assessed and one might run a risk of interpreting the unknown hadronic uncertainties as signals of physics BSM. We should also emphasize that 2HDM alone cannot explain R SM K > R exp K , and in this paper we will ignore the channels with electrons in the final state. A study along the line we are pursuing here has been initiated in Ref. [17] in which the authors computed the Wilson coefficients in the Aligned 2HDM (A2HDM), for the operators relevant to the B s → µ + µ − decay. In this paper we revisit their computation and extend it to include the operators that are needed for the phenomenological analysis of B → K ( * ) + − and other similar decays. While we broadly agree with the results of Ref. [17], there are a couple of points in which we disagree. We will examine those points, compute the remaining Wilson coefficients and use our results to discuss the phenomenological consequences on the 2HDM scenarios by comparing high−q 2 with their experimental values. We will then discuss the consequences on the similar decays with τ -leptons in the final state.
The remainder of this paper is organized as follows: In Sec. 2 we remind the reader of the main general constraints on the spectrum of scalars in 2HDM and perform a scan of parameters by assuming the lowest CP-even Higgs state to be the one measured at LHC. In Sec. 3 we write the low energy effective theory and present our results for all the Wilson coefficients in Sec. 4. We compare our results with the existing ones (in the limits in which the comparison can be made) in Sec. 5 and elucidate the subtleties related to the matching procedure in the between the full (2HDM) and effective theories in Sec. 6. Phenomenological discussion is made in Sec. 7 and Sec. 8. We briefly conclude in Sec. 9.

General constraints on 2HDM
In this Section we remind the reader of the basic ingredients of 2HDM, enumerate the parameters of the model and list the main general constraints on the spectrum of scalars which are then used to perform a scan of allowed parameters to obtain the allowed ranges of the Higgs masses and couplings.

2HDM
We consider a general CP-conserving 2HDM with a softly broken Z 2 symmetry. The most general potential can then be written as where the term proportional to m 2 12 accounts for the soft breaking of Z 2 . 2 The scalar doublets Φ a (a = 1, 2) can be parameterized as possible and they are called Type I, II, X (Lepton Specific) and Z (Flipped) 2HDM [14]. 3 To be more specific, we first write the Yukawa Lagrangian as where u and d stand for the up-and down-type quark, is a lepton flavor, f stands for a generic fermion, V for the Cabibbo-Kobayashi-Maskawa (CKM) matrix, and P L,R = (1 ∓ γ 5 )/2. A specific choice of parameters ζ f corresponds to the above mentioned types of 2HDM, which we also summarize in Table 1. Notice that the couplings ξ ϕ 0 i f appearing in the neutral Lagrangian part can be mapped onto the charged ones via Model Type X (lepton specific) cot β cot β − tan β Type Z (flipped) − tan β cot β cot β Table 1: Couplings ζ f in various types of 2HDM.

General Constraints and Scan of Parameters
To perform a thorough scan of scalars in a general 2HDM we use the general constraints summarized below.
• Stability: To ensure that the scalar potential is bounded from below, the quartic couplings should satisfy the relations [13] λ 1,2 > 0, λ 3 > −(λ 1 λ 2 ) 1/2 , and 3 The model that we call Type Z or Flipped 2HDM is sometimes referred to as Type Y.
Furthermore, the stability of the electroweak vacuum implies that which then allows us to express m 2 11 and m 2 22 in terms of the soft Z 2 breaking term m 2 12 and the quartic couplings λ 1−5 . These constraints should be combined with the necessary and sufficient condition that the minimum developed at (v 1 , v 2 ) is global [21]: • Perturbative Unitarity: An important constraint on the spectrum of scalars within 2HDM stems from the unitarity requirement of the S-wave component of the scalar scattering amplitudes. That condition implies the following inequalities [22,23] |a ± |, |b ± |, |c ± |, |f ± |, |e 1,2 |, |f 1 |, where • Electroweak Precision Tests: Finally, the additional scalars contribute to the gauge boson vacuum polarization. As a result, the electroweak precision data provide important constraint. In particular the T parameter bounds the mass splitting between m H and m H ± in the scenario in which h is identified with the SM-like Higgs, cf. Ref. [24] for example. The general expressions for the parameters S, T and U in 2HDM can be found in Ref. [25]. To derive the bounds on the scalar spectrum we consider the following values and the corresponding correlation matrix [26], The χ 2 function is then expressed as where the vector of central values and uncertainties are denoted as X = (∆S, ∆T, ∆U ) and σ = (0.11, 0.13, 0.11), while the elements of the covariance matrix are obtained via σ 2 ij ≡ σ i corr ij σ j .
As mentioned above, we identify the lightest CP-even state h with the SM-like scalar observed at the LHC with mass m h = 125.09(24) GeV [19]. To forbid the dangerous decays h → AA which could over-saturate the total width of h ( Γ SM h ), we assume that m A > m h /2. Moreover, we impose the alignment condition | cos(β − α)| ≤ 0.3, in order to ensure that the couplings of h to V = W, Z remain consistent with the values measured so far, which appear to be in good agreement with the SM predictions [27]. The abovementioned constraints are then imposed onto a set of randomly generated points in the intervals: A scan of parameters consistent with the constraints listed above favors the moderate and small values of tan β ∈ (0. 2,15]. To see that the larger values of tan β cannot be discarded it is sufficient to examine Eq. (6) in the alignment limit. For that reason, and in addition to the free scan, we perform a second scan with m H ≈ |M |, which helps us probing higher values of tan β, and we then combine results of both scans. The combined results are shown in Fig. 1 in two planes, (tan β, m H ± ) and (m A , m H ± ). From the right panel of Fig. 1 we observe that the additional scalars become mass degenerate in the decoupling region (M 2 v 2 ), as it can be easily deduced from Eqs. (6)- (9). We should also emphasize that the results of our scans agree with what has been previously reported in the literature, cf. [28].
In Sec. 8 we will confront the points allowed by our scan with the experimental measurements of exclusive b → s decays.

Effective Hamiltonian
The most general effective Hamiltonian describing the b → s transitions, made of dimension six operators, is given by [29]  where and O 7 = e/(4π) 2 m b (sσ µν P R b)F µν is the electromagnetic penguin operator. The operators with a flipped chirality, O 7,9,10,S,P , are obtained from O 7,9,10,S,P by replacing P L ↔ P R in the quark current. The dimension six operators appearing in Eq. (21) are sufficient to match the oneloop amplitude when the external fermion momenta are neglected. This, however, is not true if the computation is made with external momenta different from zero which is, in general, necessary when dealing with (pseudo-)scalar operators. For example, in order to get a correct expression for the Wilson coefficient C P one needs to consider the external momenta, which then can give rise to the contributions coming from the dimension-seven operators. One class of such terms can be related to the operators of basis (22) by equations of motion. For example, A complication arises when encountering the operators with insertion of / p b + / p s in the leptonic current, with the convention b(p b ) → s(p s ) − (p − ) + (p + ), where we also use q = p b − p s = p + + p − . A way to deal with that, adopted in Ref. [17], consists in setting p s = 0, so that / p b + / p s = / q + 2 / p s = / q = / p + + / p − , and in this way one can again, like in the previous example, use the equations of motion. That way to deal with the problem in hands, however, leads to a wrong expression for C P , for example. If, instead, one keeps all the momenta non-zero, we get a correct result. At this point we just emphasize that the matching should be performed by keeping all the external momenta different from zero and the contributions stemming from dimension-seven operators can be neglected at the very end of computation. We further elucidate this problem in Sec. 6 where we also propose a general framework for the appropriate matching between the full and effective theories in a case in which the (pseudo-)scalar bosons are explicitly taken into account.

Wilson Coefficients
After unambiguously matching the full with the effective theories we obtain the one-loop expressions for the Wilson coefficients generated by the additional scalar particles. We summarize our results in this Section. For clarity we will write them as, where the superscripts denote the types of diagrams that contributes to a given Wilson coefficient, namely, the box diagrams, the γ, Z-penguins and the (pseudo-)scalar penguins. These coefficients should be added to the (effective) ones obtained in the SM: C 7 = −0.304, C 9 = 4.211, C 10 = −4.103, and C S,P 0 [30]. 4 Henceforth, we neglect the s-quark mass and give all our results in the unitary gauge. To check the consistency of our formulas, we also performed the computation in the Feynman gauge. In the remainder of this Section we present our resulting expressions for each of the coefficients appearing in Eqs. (28)- (30). We use the standard notation, where q ∈ {b, t}, and ϕ 0 i ∈ {h, H, A}.

γ-penguins in 2HDM
The γ-penguin diagrams induced by the charged Higgs are shown in Fig. 2. The offshell and on-shell contributions can be matched onto the Wilson coefficients C 7 and C 9 , respectively, we obtain, and The dominant terms in both C NP,γ 7 and C NP,γ 9 come from the top quark contribution and are proportional to |ζ u | 2 . The terms proportional to ζ * u ζ d are suppressed by m 2 b , thus indeed subdominant.

Z-penguins in 2HDM
The Z-penguin diagrams contribute significantly to the Wilson coefficients C P , C 9 and C 10 through the diagrams shown in Fig. 3. The leading order expressions for C 9 and C 10 read, Similarly, for C P we obtain,

Charged Higgs Boxes in 2HDM
The box diagrams, peculiar for 2HDM, are drawn in Fig. 4. At low-energy they contribute to the Wilson coefficients C S and C P as, and In addition to C NP, box S,P , the tensor and (axial-)vector operators receive contributions but suppressed by the lepton mass, i.e. by x = m 2 /m 2 W . These coefficients are negligible even for decays with τ 's in the final state as it can be verified by using the expressions we provide in Appendix C.2.

Scalar penguins in 2HDM
We now turn to the effective coefficients C NP, A P , C NP, h S and C NP, H S , generated by the scalar penguin diagrams shown in Fig. 5. We recall that the total ultraviolet divergence coming from these diagrams is proportional to the factor (1 + ζ u ζ d )(ζ u − ζ d ), which vanishes due to the Z 2 symmetry (cf. Table 1). 5 The penguins with the CP-odd Higgs give rise to, where we used that ζ f ∈ R, and (1 + ζ u ζ d )(ζ u − ζ d ) = 0. Similarly, the penguins with the CP-even Higgs lead to: where λ

Comparison with Other Computations
In this Section we compare our Wilson coefficients with the results obtained in previous studies. Before doing so we should emphasize the novelties of the present work: (i) The result for C 9 in a general 2HDM with a Z 2 symmetry is new; (ii) The subleading terms O(m b ) to C 9,10 have been neglected in the previous computations, and they are included here; (iii) We provided an independent computation of the coefficients C S and C P , and elucidate inconsistencies present in Ref. [17], cf. Sec. 6 where we propose a general prescription for matching procedure when the external momenta are not neglected.
The effective coefficients C S and C P , in the context of Type II 2HDM, were first computed in Refs. [31][32][33][34][35][36]. In these papers tan β was assumed to be very large, which considerably simplifies the computation because in that case only the box diagrams give significant contributions. We agree with these results if we keep only the leading terms in tan β in our expressions, namely, Along the same lines, the leading order QCD corrections to the same coefficients were included in Ref. [37]. Recently, the computation of C S and C P was extended to the context of a general A2HDM, which comprises all four types of 2HDM with Z 2 symmetry discussed  here but without the usual assumption of large tan β [17]. We agree with their general results, except for the expression for C NP, Z P which differs from the one reported in the present paper. The disagreement comes from the fact that the authors of Ref. [17] worked with the assumption p s = 0, which appears not to be fully appropriate. 6 By keeping p s = 0 one realizes that the computation of Z-penguin leads to two independent terms, one proportional to p H = p b +p s and the other to q = p b −p s . By using equations of motion, C P,S correctly receive contributions from the terms proportional to q, but not from those proportional to p H . With p s = 0 only one invariant appears, because p H ≡ q, and thus the resulting C P,S also receive spurious contributions from p H .
Regarding the other Wilson coefficients, the first computations of C 7 for a general 2HDM have been performed in Ref. [38], then in Refs. [39,40] and [41] where the leading and subleading QCD corrections were included too. Our results are consistent with those, as well as with the expression for C 10 presented in Ref. [35] and more recently in Ref. [17]. The only difference with respect to those results is that we include the subleading terms in m b .

Matching Procedure
In this section we discuss in more detail the matching of the one-loop amplitudes when the nonzero external momenta are considered. We stress once again that keeping external momenta non-zero is necessary to obtain the correct values for the Wilson coefficients C S,P . As we mentioned in Sec. 3 the insertion of external momenta result in dimensionseven operators which can be simplified by using equations of motion, except in the cases when the lepton momenta are to be contracted with the quark current and/or the quark momenta to be contracted with the lepton current. The amplitudes which need a special treatment, to leading order in external momenta, are: where i, j = L, R and s, b, are the fermion spinors. Note again that our convention is b(p b ) → s(p s ) − (p − ) + (p + ), and q = p b − p s = p + + p − . In order to keep our discussion general, we first extend the Hamiltonian (21) and include the following operators 6 We should emphasize that we were able to reproduce the expression for C NP, Z P reported in Ref. [17] by taking p s = 0, which however is not an appropriate assumption as we argue in the text. where with i, j = L, R. 7 We reiterate that even though these operators are suppressed by 1/m W , they are necessary to unambiguously match the loop induced amplitudes with the effective field theory. The above choice of the basis of dimension-seven operators is convenient since they do not contribute to B(B s → µ + µ − ), while for the other decays their hadronic matrix elements are easy to calculate. By using the Fierz rearrangement and by applying the field equations, the amplitudes (42) are reduced to To remain completely general, in the above equations we also kept the lepton mass and the mass of s-quark different from zero. As an example we show the validity of Eq. (47). Using p − − p + = 2p − − q, and by the multiple use of field equations, we can write: By applying the Fierz identity once again, we arrive at, Clearly, for the appropriate matching of these amplitudes to the effective theory, the operators appearing in Eq. (21) are not enough and the extended basis given in Eq. (43) is necessary. Once the matching is performed, the operators from Eq. (21) could be neglected since they are 1/m W suppressed with respect to the dominant (dimension six) ones. This delicate point can then be verified explicitly by computing the Wilson coefficients C T q RL and C T q RR which come from the Z-penguin diagrams and the coefficients C T LL = (C T LR ) * generated by the box diagrams. Their explicit expression is given in Appendix C.1.
We can now easily understand the source of our disagreement with Ref. [17]. If one sets p s = 0 in A q RR of Eq. (42), then just like in Ref. [17] one could write / p b + / p s = / p b = / q which, by means of equations of motion, yields which then in the actual computation gives a contribution to C P . With our procedure, we understand that this contribution does not come from C P but actually from √ x C T q RL . In other words, and by using our definition of operators and of the effective Hamiltonian, we Therefore the Wilson coefficient C P of Ref. [17] contains the Wilson coefficient of the operator O T q RR , the matrix element of which is not equal to the matrix element of the operator O P but is, instead, suppressed by m W as we explicitly check in the next section. For that reason the Wilson coefficient of Ref. [17] is not correct.
In this Section we give the expressions for B(B s → µ + µ − ) and B(B → Kµ + µ − ) to which we also include the contributions of the operators given in Eq. (44). Those additional operators were necessary for the appropriate matching procedure between the full and the effective theories. However, since they are suppressed by 1/m W they are expected to be negligible with respect to the dominant operators entering the effective Hamiltonian (21). The purpose of this exercise is to check whether or not the size of the matrix elements of the operators (44) is indeed numerically insignificant for phenomenology.

B s → µ + µ −
On the basis of Lorentz invariance and invariance of the strong interaction with respect to parity, one can easily verify that B s → µ + µ − is not affected by the operators O T q i,j and O T i,j , with i, j = L, R. The expression for the decay rate of this process remains the standard one where β = 1 − 4m 2 /m 2 Bs . To compare Eq. (57) with the available experimental value, one needs to take into account the effects of B s − B s oscillations which, to a good approximation, amounts to [42] where y s = ∆Γ Bs /(2Γ Bs ) = 0.061 (9), experimentally established by the LHCb Collaboration [43]. As we mentioned before, the dimension-seven operators (44) were chosen in such a way that they do not contribute the B s → + − decay amplitude.

B → Kµ + µ −
In contrast to B s → + − , the decay B → K + − receives contributions from the operators of the extended basis (44). To write the decay amplitude in a compact form, it is convenient to use the formalism of helicity amplitudes (HA's). In the absence of the (pseudo-)scalar operators, the total amplitude can be schematically written as By describing the decay mode as B → KV * → K + − , where V * is a virtual vector boson, one can decompose the total decay amplitude in terms of HA's, where ε µ V (m) (with m, n = 0, t, ±) are the V * -boson polarization vectors, explicitly defined in Appendix A. We repeat that the above decomposition is valid as long as the scalar and the pseudoscalar operators are not present. To incorporate those contributions unambiguously one can assume the lepton masses to be unequal (m 1 = m 2 ) and then apply the Ward identities, to absorb the (pseudo-)scalar terms in the time-like coefficients A L(R) t . By taking the limit m 1 = m 2 in the final expression one ends up with the desired HA's and the total decay amplitude, from which is then easy to compute the decay rate [44]. Notice that the contributions from C ( ) S,P enter the amplitudes A S and A t defined as, More details regarding this point can be found in Ref. [44]. We also need to stress that all the helicity amplitudes are the q 2 -dependent functions, A i ≡ A i (q 2 ). By applying the method briefly sketched above we obtain, where the explicit expressions for the helicity amplitudes are: where the normalization factor also accounts for the remaining phase space, namely, For shortness, in the above formulas, we used λ q = λ( q 2 , m , m ) and The kinematic conventions and the form factor definitions are collected in Appendix A. In the limit in which the derivative operators vanish we retrieve the usual expression for differential branching fraction [44]. The choice of dimension-seven operators (44) is convenient also because their matrix elements are proportional to the original hadronic matrix elements multiplied by iq µ . As it can be seen from the above expressions the coefficients C T i,j and C T q i,j enters the above formulas with the explicit 1/m W -suppression factor. In other words, with the above formulas and by using the Wilson coefficients presented in the previous Sections, we see that the derivative operators (44) are indeed irrelevant for phenomenology. Their presence is therefore essential for the unambiguous matching procedure in the computation of Wilson coefficients but they do not alter the phenomenological analysis even at the sub-percent level.

Phenomenology and discussion
In this Section we use our results for Wilson coefficients and compare the experimental data for the exclusive b → s + − modes with various types of 2HDM. We decided to focus on B(B s → µ + µ − ) exp = (2.8 +0.7 −0.6 ) × 10 −9 [1], and B(B → Kµ + µ − ) exp high q 2 = (8.5 ± 0.3 ± 0.4) × 10 −8 [3], where "high q 2 " means that the decay rate has been integrated over the interval q 2 ∈ [15, 22] GeV 2 . The reason for opting for these decay modes is that the relevant hadronic uncertainties are under good theoretical control. The hadronic quantity entering the B s → µ + µ − decay amplitude is the decay constant, f Bs . It has been abundantly computed by means of numerical simulations of QCD on the lattice (LQCD) and its value is nowadays one of the most accurately computed hadronic quantities as far as B (s) -mesons are concerned [16]. The hadronic form factors entering the B → Kµ + µ − decay amplitude have been directly computed in LQCD only in the region of large q 2 's [46,47], which explains why we use B(B → Kµ + µ − ) exp high q 2 to do phenomenology. Furthermore, since the bin corresponding to q 2 ∈ [15, 22] GeV 2 is rather wide and away from the very narrow charmonium resonances, the assumption of quark-hadron duality is likely to be valid [45].
By using the recent LQCD results for the form factors provided by HPQCD [46] and MILC Collaborations [47], the SM results are both being about 2σ larger than the experimental value measured at LHCb. 9 Since the current disagreement between theory and experiment needs to be corroborated by more data, we decided to impose all the constraints to 3σ accuracy. We will then discuss the impact of B(B → Kµ + µ − ) exp high q 2 on 2HDM if the current discrepancy remains, i.e. by requiring the 2HDM to compensate the disagreement between theory (SM) and experiment at the level of 2σ and more. Notice also that the measured B(B s → µ + µ − ) exp is slightly smaller than predicted, B(B s → µ + µ − ) SM = (3.65 ± 0.23) × 10 −9 [2].
We now use the results of our scan from Sec. 2.2, require the 3σ agreement between experiment and theory, which means that we add the generic 2HDM Wilson coefficients derived in the previous Section to the SM values. The result, in the plane (tan β, m H ± ), is shown in Fig. 6 for each type of 2HDM discussed in Sec. 2. We learn that both B(B s → µ + µ − ) and B(B → Kµ + µ − ) high q 2 exclude the low tan β 1 region regardless of the type of 2HDM considered. The limit of exclusion of low tan β coming from B(B → Kµ + µ − ) high q 2 is slightly larger than the one arising from B(B s → µ + µ − ). The limit on low tan β obtained in this way for each of our four models is given in Tab. 2.
Besides excluding tan β 1, it may appear as a surprise that the large tan β are not excluded by these data. The reason for that is the fact that the (pseudo-)scalar Wilson coefficient, with respect to the dominant (axial-)vector one, comes with a term proportional to (m Bs /m W ) 2 which suppresses the large tan β values. This feature can be easily verified in the Type II model for which the coefficients C S,P , in the large tan β limit, are given in Eq. (41). This is why only a small number of points have been eliminated from our scan of Type II model at large tan β but relatively light m H ± .
Model Type I Type II Type X Type Z tan β > 1.0 > 0.9 > 1.0 > 0.9 Table 2: Allowed values of low tan β (at 99% CL) for the different 2HDMs. See text for details.
Since the SM value is in slight tension with B(B → Kµ + µ − ) exp high q 2 at the 2.1σ level, we can now check which of the models discussed in this paper can be made consistent with the experimental data if any disagreement beyond 2σ between theory (SM) and experiment is to be attributed to 2HDM. It turns out that two such models are Type II and Type Z 2HDM, which we illustrate in Fig. 7. For the other two scenarios (Type I and Type X) the NP contributions are either too small or already in conflict with B(B s → µ + µ − ) exp . From Figs. 7 and 8 we see that in order to explain the discrepancy one needs a relatively light charged scalar: (i) m H ± 735 GeV and tan β > 2.3 in the Type II scenario, and (ii) m H ± 380 GeV and tan β > 3.5 for the Type Z scenario. Since the masses of the additional scalars are correlated, we see that m H and m A become bounded as well, cf. Fig. 8. In the case of Type II and Type Z 2HDM an additional bound on the charged Higgs has been recently derived from the inclusive mode B(B → X s γ). After comparing the experimental spectra with theoretical expressions in which the higher order QCD corrections have been included, the lower bound m H ± > 570 GeV (95% CL) was obtained in Ref. [48] (c.f. also Ref. [49]). This bound is superposed on our results in Figs. 7 and 8, which then also eliminates Type Z 2HDM. Furthermore, we can say that the requirement of agreement between theory and experiment to 2σ, for the quantities discussed in this Section, reduces the available space of parameters for Type II 2HDM to m H ± ∈ (570, 735) GeV, and tan β ∈ (16,35), while the available range of values for the mass of the CP-odd Higgs becomes m A ∈ (145, 865) GeV. Figure 7: Results of the scan in Fig. 1 after imposing the b → s constraints to 2σ accuracy. The hatched area is excluded by B(B → X s γ) at 95% [48]. See Fig. 6 for the color code. In what follows we will assume that the 2σ disagreement of the measure B(B → Kµ + µ − ) exp high q 2 with respect to the SM prediction indeed remains as such in the future and discuss the consequences on the decays B(B s → τ + τ − ) and B(B → Kτ + τ − ) high q 2 if the Type II 2HDM is used to explain the disagreement. From Eq. (57) we can see that where the only remaining m dependence comes from the last numerator in the factor multiplying |C S − C S | 2 in Eq. (57). In Fig. 9 we illustrate the validity of the above equality. Notice that a tiny departure from equality comes from the large tan β values which enhance the C S contribution. In other words, the current experimental result B(B s → µ + µ − ) exp , which is slightly lower than the one predicted in the SM, is expected to lead to B(B s → τ + τ − ) exp compatible or slightly lower than predicted in the SM. The cancellation of the lepton mass in B(B s → + − ) 2HDM , discussed above, does not occur in B(B → K + − ) 2HDM high−q 2 . As a result we obtain, where we omitted the subscript "high-q 2 " to avoid too heavy a notation. Illustration is provided in Fig. 9. We can rephrase this observation with an equivalent statement: To be fully explicit, we obtain

Conclusion
In this paper we computed the leading order Wilson coefficients relevant to the exclusive b → s + − decays in the framework of 2HDM with a softly broken Z 2 symmetry. Most of these Wilson coefficients have been computed previously but in the limit of large tan β, which we extend here to a generic setup. We also included O(m b ) corrections, which were neglected in the previous computations. Regarding the (pseudo-)scalar Wilson coefficients, we elucidated the issue of unambiguous matching of the one-loop amplitudes between the full and effective theories which requires extending the basis of operators in the effective theory by including two types of operators suppressed by 1/m W (altogether, eight new operators). We pointed out that for the appropriate identification of the Z-penguin contribution to the Wilson coefficient C P it is necessary to keep all external momenta different from zero. After having computed the full set of Wilson coefficients we were able to make a phenomenological analysis by focusing on B(B s → µ + µ − ) and B(B → Kµ + µ − ) high−q 2 , the quantities which are measured at LHC and for which the hadronic uncertainties are under good theoretical control (computed in LQCD). After carefully scanning the parameter space of 2HDM with a softly broken Z 2 symmetry, we tested various types of 2HDM against the experimental data for B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp high−q 2 , and found that to 3σ the values of low tan β 1 are excluded for all types of 2HDM's considered here.
If, instead, we require the 2σ agreement with experiment, then only Type II and Type Z models provide a viable description of the data. After combining ours with the bound on the charged Higgs deduced from the inclusive b → sγ decay, we find that the Type Z model can be discarded and Type II : m H ± ∈ (570, 735) GeV, m A ∈ (145, 865) GeV, tan β ∈ (16,35). (77) We also discussed the repercussions of the current results on the decays B(B s → τ + τ − ) and B(B → Kτ + τ − ) high−q 2 .
In the dilepton rest frame the components of the leptonic four-vectors are given by p µ − = (E , |p | sin θ , 0, |p | cos θ ), p µ + = (E , −|p | sin θ , 0, −|p | cos θ ), where E = q 2 /2, and θ is the angle between − (in the dilepton rest frame) and the line of flight of the two leptons (in the B-meson rest frame). The momentum p is given by

Hadronic matrix elements
For completeness we also give the definitions of the decay constant (f Bs ) and of the form factors [f +,0,T (q 2 )], quantities which parametrize the hadronic matrix elements relevant to the processes discussed in this paper: where for B → K + − the kinematically accessible q 2 values lie in the interval 4m 2 ≤ q 2 ≤ (m B − m K ) 2 . Notice that we do not write explicitly the scale dependence of the quark masses, nor of the scalar and tensor densities and of the form factor f T (q 2 ). In the actual computations the MS values of these quantities are taken at µ = m b .
charges of the initial particles. For the trilinear scalar interactions, we have where the trilinear couplings read These results agree with the ones given in Refs. [17,51] after the appropriate change of basis and/or conventions. 10

C Scalar penguins and auxiliary functions
In this Appendix we give the expressions for the Wilson coefficients generated by each diagram shown in Fig. 5. We also give the expressions for the auxiliary functions (f i and g i ) used in this paper. The penguins arising from coupling to ϕ 0 i ∈ {h, H, A} contribute to the effective coefficient C S,P and can be generically written as where C k,ϕ 0 i is a common coefficient generated by the diagram k, with k = 1, . . . , 18. Since, in our framework, ζ h , ζ H ∈ R and ζ A ∈ i × R, it is clear that the CP-even scalars h and H contribute only to C S , whereas the CP-odd Higgs A contributes only to C P , as expected from the assumption of CP conservation. We obtain in the unitary gauge, + these contributions here. We obtain: and C NP, box .