Revisiting $B_s \to \mu^+\mu^-$ in the two-Higgs doublet models with $Z_2$ symmetry

We revisit the rare leptonic decay $B_s \to \mu^+ \mu^-$ in the two-Higgs doublet models with a softly broken $Z_2$ symmetry, namely type-I, type-II, type-X and type-Y 2HDMs. We have derived the relevant full one-loop Wilson coefficients of the four 2HDMs from the recent calculation in the aligned two-Higgs doublet model by Li, Lu and Pich, which could be mapped to all the four 2HDMs for both large and small $\tan\beta$. It is found that, a new term associated with the soft $Z_2$ symmetry breaking parameter $M$ can be enhanced by $\tan^2\beta$ in the type-II 2HDM, which has not been considered in the literature. Imposing both theoretical and experimental constraints, we have renewed the bounds on the parameter spaces of the four 2HDMs. Different from our previous paper, however, we find that all the four 2HDMs give sizable and similar contributions to $\overline{\mathcal B}(B_s \to \mu^+ \mu^-)$ within the stringently restricted parameter spaces, but very tiny as regards the mass-eigenstate rate asymmetry $\mathcal A_{\Delta\Gamma}$; this makes it unfeasible to discriminate the four types of 2HDM with the correlations between the observables in $B_s \to \mu^+ \mu^-$ decay.


Introduction
The discovery [1,2] of a new boson with a mass close to 125 GeV has been well anticipated as the standard model Higgs boson [3][4][5] and provided the first experimental evidence of the Higgs mechanism [6][7][8]. It is a great triumph, but not an end, of the giant campaign for Higgs hunting in the development of particle physics. Although the subsequent more precise measurements [9][10][11][12][13] at the LHC have shown the properties of the Higgs boson are well consistent with the predictions of the standard model (SM), the precision of the current experimental data still leave open the possibility of an extended Higgs sector [14,15]. Among many new physics scenarios beyond the SM, the two Higgs doublet models (2HDM) [16][17][18] are the simplest extensions of the SM.
Among the rare B-meson decays, the leptonic processes B q → µ + µ − (q = d, or s) are of special interest [49,50]. They suffer from very few hadronic uncertainties and are induced by FCNC transitions, which make them sensitive probes to the effects of physics beyond the SM, especially models with a non-standard Higgs sector [51][52][53][54][55]. Recently, the next-to-leading order (NLO) electroweak corrections and the next-to-next-to-leading order (NNLO) QCD corrections [56][57][58] in the SM have been calculated. On the BSM side, a full one-loop calculation in the aligned 2HDM (A2HDM) has been performed in ref. [59].
Motivated by this progress, in this paper we perform a detailed study of the B s → µ + µ − decay within the 2HDMs with Z 2 symmetry. At present, this process is calculated in the type-II 2HDM in large tan β limit only [60][61][62]. Using the Higgs base correspondence between the A2HDM and the 2HDMs, we will derive the relevant full one-loop Wilson coefficients of the four variant 2HDMs contributing to the B s → µ + µ − decay from the recent A2HDM results [59] without the large tan β approximation. We also investigate the possibility to discriminate the four different types of 2HDM in the light of the recent collider and flavor physics data, as an update of our previous work [63]. We combine the constraints from B s,d → µ + µ − , B s,d −B s,d mixing, B → τ ν andB → X s γ [64,65], with the experimental data from the direct search for Higgs bosons at LEP [66], Tevatron [67,68] and LHC [69,70], and the constraints from perturbativity, tree-level vacuum stability and perturbative unitary. For the B s → µ + µ − decay, the correlations between its branching ratio and the mass-eigenstate rate asymmetry A ∆Γ are also reevaluated with the constrained parameter space of the 2HDMs obtained in this paper.
We have found that A ∆Γ can slightly deviate from the SM prediction in the type-II 2HDM only, and that the ratio of time-integrated B(B s → µ + µ − ) gets similar contributions from the four 2HDMs; this makes it very hard to discriminate the four types of 2HDMs with the correlation between A ∆Γ and B(B s → µ + µ − ) as suggested in our previous work [63].
The paper is organized as follows. In section 2, we give a brief overview of the B s → µ + µ − decay. In section 3, full one-loop contributions from the 2HDMs with Z 2 symmetry are derived explicitly. In section 4, we give our detailed numerical results and discussions. We conclude in section 5. The relevant theoretical formulas are recapitulated in the Appendix.
2 B s → µ + µ − in the SM In the SM, the leptonic decays B q → µ + µ − (q = d or s) arise from the W box and Z penguin diagrams. Generally, these decays can be described by the low-energy effective Hamiltonian where α e denotes the QED fine-structure constant and V ij the CKM matrix elements. The semi-leptonic operators are defined as In the SM, the contributions from the scalar operators O S and O P are highly suppressed (the corresponding Wilson coefficients are given in eq. (A.1).), but C 10 will play the dominant role.
This progress will be incorporated into our calculations.
With the effective Hamiltonian eq. (2.1), the branching ratio of B q → µ + µ − reads where m Bq , τ Bq and f Bq denote the mass, mean lifetime and decay constant of B q meson respectively. The short-distance contributions S and P are defined as As discussed in the following section, there is no BSM phase in the 2HDMs with Z 2 symmetry.
Therefore, we only consider the case that both S and P are real in this paper.
As pointed out in ref. [74], the measured branching ratio of B q → µ + µ − should be the time-integrated one, denoted by B(B q → µ + µ − ). In order to compare with the experimental measurements, the sizable effect of B s −B s oscillations should be taken into account [74,75], and one has where the mass-eigenstate rate asymmetry A ∆Γ can be expressed as The observable A ∆Γ is independent of the branching ratio of B s → µ + µ − and provides complementary information on the short-distance structure of this decay. In the SM, A ∆Γ = +1.
Following ref. [74], it is convenient to introduce the ratio where both hadronic uncertainties and CKM matrix elements are canceled out.
3 B s → µ + µ − in the 2HDMs with Z 2 symmetry In the 2HDMs with Z 2 symmetry, b → sµ + µ − processes receive contributions from box diagrams with charged Higgs and penguin diagrams with Z boson and neutral Higgs bosons. The Wilson coefficient C 10 has been calculated in the type-II 2HDM [53]. For C S and C P , only the leading contributions in the large tan β limit have been computed in the type-II model [60][61][62]. However, the remaining contributions could be important for some specific tan β values in the other types of 2HDMs. In this section, we first of all give a brief introduction to the 2HDMs with

2HDMs with Z 2 symmetry
The 2HDM extends the SM Higgs sector with an additional scalar doublet. With the two Higgs doublets Φ 1 and Φ 2 , the CP-conversing 2HDM potential with a softly broken Z 2 symmetry reads [18] is a soft Z 2 symmetry breaking term and the parameters m 1−3 and λ 1−5 are real. The two Higgs doublets Φ 1 and Φ 2 can be generally parameterized as where the two vacuum expectation values (vev) v 1 and v 2 are real and positive. From the vacuum condition [76] they can be expressed as other parameters in the Higgs potential, where λ 345 = λ 3 + λ 4 + λ 5 is defined. By introducing the vev v (v = v SM = 246 GeV), the mixing angle β and the soft Z 2 symmetry breaking parameter M as v 1 = v cos β, v 2 = v sin β and M 2 = m 2 3 /s β c β , we can use (v, β, M, λ 1−5 ) as independent 2HDM potential parameters.
Physical Higgs states are obtained by the following rotations: where the rotation matrix is given by The mixing angle α is determined by the Higgs potential of eq. (3.1) [76], In the 2HDM with Z 2 symmetry, the physical Higgs spectrum consists of five degrees of freedom: two charged scalars H ± , two CP-even neutral scalars h and H, and one CP-odd neutral scalar A. The quartic couplings λ i in the Higgs potential can be expressed in terms of their masses as [76] Therefore, the eight parameters in the Higgs potential m 1−3 and λ 1−5 can be rewritten equivalently by the four physical Higgs masses m h , m H , m A , m H ± , the two mixing angles α and β, the vev v = v SM , and the Z 2 symmetry breaking parameter M . In the case of λ 1 = λ 2 , which is considered in ref [61,62], M can be eliminated and the 2HDM potential parameters can be In the interaction basis, the general Yukawa Lagrangian of the 2HDM can be written as (3.9) whereΦ i = iσ 2 Φ * i , Q L and L L denote the SM quark and lepton doublets, and u R , d R , and e R are the right-handed up-type quark, down-type quark and lepton singlet, respectively. The Yukawa coupling matrices Y u,d, i are 3 × 3 complex matrices in flavor space.
In order to avoid tree-level FCNC, a discrete Z 2 symmetry is introduced [48]. All the possible nontrivial Z 2 charge assignments are listed in table 1, which define the four well-known types of 2HDM, i.e. type-I, type-II, type-X and type-Y. In the mass-eigenstate basis, the Yukawa interactions can be written in the form where P L,R = (1 ∓ γ 5 )/2. The Yukawa couplings ξ f h,H,A in the four types of 2HDM are listed in table 2. In addition, the couplings of the light CP-even Higgs boson h to gauge bosons or ZZ can be written as g hV V = sin(β − α)g SM hV V , which is normalized to the corresponding couplings of the SM Higgs boson g SM hV V [18]. Recently, the LHC Run I data confirm the SM Higgs-like nature of the 125 GeV boson discovered at the LHC [3][4][5]. If the light CP-even Higgs h in the 2HDM is identified with the observed 125 GeV boson, global fits to the LHC Higgs data suggest that all four types of 2HDM should lie close to the so-called alignment limit [77][78][79][80][81][82][83] sin(β − α) = 1, For recent studies on the alignment limit in the 2HDM, we refer to ref. [84,85].
Since the 2HDMs with Z 2 symmetry are particular cases of the A2HDM [47], there exists a one-to-one correspondence for Yukawa couplings between these two models. However, the correspondence is not so straightforward for Higgs cubic couplings. Unlike the 2HDMs with Z 2 symmetry, the A2HDM potential is usually defined in the so-called "Higgs basis" [87], in which only one Higgs doublet gets a nonzero vev. Therefore, the parameter tan β defined in the NFC 2HDMs is not a physical parameter in the A2HDM [88].
In both the A2HDM and the NFC 2HDMs, B s → µ + µ − decay is induced by gauge boson Z, Goldstone boson G 0 , and Higgs bosons ϕ ≡ {h, H, A} penguin diagrams, as well as box diagrams mediated with W ± , H ± , and G ± . To one-loop level, their contributions to the Wilson coefficients are divided into the following different parts: . (3.14) For self-contained, we present the Wilson coefficients after the correspondences made in appendix A.  .4), and the Higgs cubic couplings are defined as where the soft Z 2 symmetry breaking parameter M has been defined in sec 3.1.
In the literature [60][61][62], it is found that the Wilson coefficients can receive large tan β The soft Z 2 symmetry breaking parameter M is associated with the spontaneous CP breaking [16,[89][90][91] and characterizes the masses of all the Higgs bosons [76]. This parameter enters the B s → µ + µ − decays through the Higgs penguin diagrams. However, it is found that the M term can not make more significant contributions than other terms of the Wilson coefficient C S . Here, we would choose h as the Higgs boson discovered by ATLAS [1] and CMS [2] and take the alignment limit β − α = π/2, which is favored by the current 2HDM fits [77][78][79][80][81][82][83]. Then the cubic couplings in eq. (3.16) read (3.17) Focusing on the coupling λ h H + H − , it can be seen from eqs. (2.4) and (3.17) that large contributions from this coupling would require |M 2 − m 2 H ± |/v 2 m 2 W /m 2 B . However, we know |M 2 − m 2 H ± |/v 2 = |λ 4 + λ 5 |/2 < 4π from the 2HDM vacuum condition [76] and perturbativity [92]. It is also noted that, the Higgs penguin diagrams can be enhanced by very large tan β or cot β. In all the four types of 2HDMs, the λ h H + H − contributions could be enhanced by large cot 2 β. In practice, cot β 3 has been excluded by the perturbativity [92]. Similarly, the Among the four models, this contribution is enhanced by tan 2 β only in type-II 2HDM. However, the ratio M 2 /m 2 H still suffers from the theoretical constraints, which will be discussed with numerical results in the following section. In the particular case of the type-II 2HDM, our result of C 10 agrees with the one calculated in ref. [53]. diagrams, the calculations are performed again in the large tan β limit but with the assumption λ 1 = λ 2 for the couplings in the Higgs potential 3 . Considering only terms proportional to tan 2 β, our result agrees with the one of ref. [60] in the case of λ h H + H − = λ H H + H − = 0, and those of ref. [61,62] in the case of λ 1 = λ 2 . Generally, the 2HDM contains eight free parameters, i.e., m 1−3 and λ 1−5 in the Higgs potential of eq. (3.1). They can be rewritten equivalently in terms of the Higgs masses m h , m H , m A , m H ± , the mixing angles α and β, the parameter M , and the vev v = v SM . If the condition λ 1 = λ 2 is assumed, M can be expressed by the other parameters, as shown in eq. (3.8). It is the reason why terms depending on the Z 2 symmetry breaking parameter M were absent in the previous calculations [60][61][62], but are present in this paper.

Numerical Analysis
Searches for B s,d → µ + µ − decays have been performed at the BaBar, Belle, and Tevatron (for a review, see ref. [97]). At the LHC, measurements by CMS [99] and LHCb [100]  Thus, strong constraints on the 2HDM parameters are expected.
In the NFC 2HDMs, the relevant parameters are the two mixing angles α and β, four Higgs mass parameters m H ± , m h , m H , and m A . In the B s,d → µ + µ − decays, the Z 2 symmetry agree with each other after the erratum for ref. [61] has been taken into account. In addition, there is a typo in eqs. (3.30) and (3.31) of ref. [61]: a global factor α e /π should be included. 3 In ref. [61,62], the convention for the Higgs potential (i.e., the couplings λ i ) is different from the one defined in eq. (3.1). This condition is also expressed as λ 1 = λ 2 by the couplings used in our paper.
breaking parameter M also enters into the decay amplitude and is independent from these parameters. As discussed in ref. [102,103], we choose the light neutral Higgs h in the 2HDM as the SM Higgs observed at the LHC and adopt the alignment limit sin(β − α) = 1. Then the model parameters are reduced to (m H , m A , m H ± , M, tan β). As discussed in ref. [63], we shall restrict these parameters in the following ranges: Starting from these parameter spaces, we will start our numerical scan.
In the numerical analysis, we impose experimental constraints in the same way as in ref. [63].
After considering the current experimental data, the allowed parameter spaces of all the four 2HDMs are obtained. Since the constraints from B d → µ + µ − appear to be more or less weaker than those from B s → µ + µ − , we only show the results from the latter one, which are plotted in the (tan β, m H ± ) plane in figure 1(b). Compared to our previous results [63], the parameter space with small tan β is excluded for all the four types of 2HDMs. This change is Combining all the constraints aforementioned, we obtain the survived parameter space of all the four types 2HDMs, as an update of our previous results [63], which is shown in the (tan β, m H ± ) plane in figure 2(a). It is found that the small tan β region is restricted for all the four models by B s −B s mixing and B → X s γ , while the large tan β region is constrained only in the type-II 2HDM by B → τ ν and B s → µ + µ − decays. Compared to our previous results, the current constraints on the large tan β region in the type-II 2HDM are more stringent. This is mainly because the theoretical constraints are included in the current analysis.
In these constrained parameter spaces of the four 2HDMs, the correlations between the observables A ∆Γ and R defined in eqs. (2.6) and (2.7) are reevaluated, which are presented in figure 2(b). Unlike our previous results [63], the correlations in the four different types of 2HDMs are almost indistinguishable. The allowed ranges of R are the same for all the four models, while A ∆Γ can deviate slightly from the SM prediction only in the type-II 2HDM.
It is found that the difference from our previous results is mainly caused by the theoretical constraints and the new full one-loop Wilson coefficients considered in the current analysis. In the type-II 2HDM, the bounds on tan β are more stringent compared to our previous results as discussed above. Thus, the allowed range of C S is restricted more stringently in the current analysis. In this case, A ∆Γ can deviate from the SM prediction very tiny, which can be seen from eq. (2.6). As discussed in sec. 3, our results of the Wilson coefficients can also be applied to the small tan β region in all the four models, while some terms are not included in C P used in our previous analysis. In the case of small m A , C P is enhanced and these terms make the allowed regions of R in the type-I and type-Y 2HDMs as large as the one in the type-X model.
Meanwhile, the value of R is almost independent of A ∆Γ in the type-II 2HDM.

Conclusion
In this paper, we have performed an updated analysis of the rare leptonic decay B s → µ + µ − in the 2HDM with a softly broken Z 2 symmetry. We have derived the full one-loop Wilson coefficients C 10 , C S and C P from the recent A2HDM results [59], which can be applied to the contributions of all the four types of 2HDMs for both large and small tan β value. Our main conclusions are summarized as follows: • Compared to C 10 , the Wilson coefficients C S and C P are negligible in the entire 2HDM parameter space, except for large tan β in the type-II 2HDM or small CP-odd Higgs mass m A in the four models. In addition, only the Wilson coefficients C S and C P in the type-II 2HDM can be enhanced by large tan β.
• The soft Z 2 symmetry breaking parameter M enters into the Higgs penguin diagrams and affects the Wilson coefficient C S . The dominant contributions are proportional to M 2 /m 2 H and enhanced by tan 2 β in the type-II 2HDM, which have not been considered in the literature [60][61][62]. However, after combing the theoretical constraints from perturbativity, vacuum stability and perturbative unitarity, we have found that the parameter M can not make more significant contributions than other terms in the Wilson coefficients.
• After imposing the experimental constraints, regions with small tan β are excluded for all the four types of 2HDM, which are quite different from our previous results [63]. As expected, large tan β region is only excluded in the type-II 2HDM.
As an update of our previous analysis [63], we have also investigated the possibility to

A The Wilson coefficients in the SM and the 2HDMs
In this appendix, we recapitulate the relevant expressions of the Wilson coefficients in the SM and the four types of the 2HDMs for completeness, which are obtained from ref. [59].
In the SM, the one-loop Wilson coefficients of the scalar operators can be written as In the unitary gauge, their expressions read where C h, SM S, Unitary denotes the contributions from the SM Higgs penguin diagrams. The other Wilson coefficients C box, SM S,P, Unitary and C Z, SM P, Unitary are same in the SM and the 2HDMs. In the four types of the 2HDMs, the various contributions in the Wilson coefficients of eq. (3.13) are obtained by the replacement of the Yukawa couplings in eq. (3.14), which are given in the unitary gauge, The Higgs penguin contributions C ϕ, 2HDM S,P, Unitary have been given in eq. (3.15), where the functions g (a) 0−3 are defined as Here the one-loop functions f i are abbreviated as f i ≡ f i (x t , x H ± ) with the definitions f 1 (x, y) = 1 2(y − x) x − y + y ln y − x ln x , (A.5) It is noted that the divergence in the Higgs penguin diagrams at one-loop level is canceled by a FCNC local operator in the A2HDM [59]. In the 2HDMs with Z 2 symmetry, we find that the divergence automatically vanishes after adding all the Higgs penguin contributions.