$B_{s,d}^0 \to \ell^+\ell^-$ Decays in Two-Higgs Doublet Models

We study the rare leptonic decays $B_{s,d}^0 \to \ell^+\ell^-$ within the general framework of the aligned two-Higgs doublet model. A complete one-loop calculation of the relevant short-distance Wilson coefficients is presented, with a detailed technical summary of the results. The phenomenological constraints imposed by present data on the model parameters are also investigated.


Introduction
The discovery [2,3] of a Higgs-like boson at the LHC has placed the last missing piece of the Standard Model (SM), which is one of the greatest achievements of modern particle physics. However, it is widely believed that the SM cannot be the fundamental theory up to the Plank scale, and many theories beyond the SM (BSM) claim that new physics (NP) should appear around the TeV scale.
One of the simplest extensions of the SM is the addition of an extra Higgs doublet [4]. Two scalar doublets are present in several BSM theories, for instance in supersymmetry. Two-Higgs doublet models (2HDMs) with generic Yukawa couplings give rise to dangerous tree-level flavour-changing neutral currents (FCNCs) [5]. This can be avoided imposing discrete Z 2 symmetries [6] or, more generally, assuming the alignment in flavour space of the two Yukawa matrices for each type of right-handed fermions [7]. 1

Speaker
The leptonic decays B 0 s,d → + − play a very special role in testing the SM and probing BSM physics. They are very sensitive to the mechanism of quark flavour mixing, and their branching ratios are extremely small due to the loop suppression and the helicity suppression factor m /m b . Since the final state involves only leptons, the SM theoretical predictions are very clean [8]: which include next-to-leading order (NLO) electroweak corrections [9] and next-to-next-to-leading order (NNLO) QCD corrections [10]. The weighted world averages of the CMS [11] and LHCb [12]   It is convenient to define the 2HDM in the "Higgs basis" where only one scalar doublet gets a nonzero vacuum expectation value v = ( √ 2G F ) −1/2 246 GeV: The first doublet contains the Goldstone fields G ± and G 0 . The five physical degrees of freedom are given by the two charged fields H ± (x) and three neutral scalars The latter are related with the S i fields through an orthogonal transformation R, which defines the neutral mass eigenstates: The mass matrix M of the neutral scalars is fixed by the scalar potential: , where µ 1 , µ 2 and λ 1,2,3,4 are real, while µ 3 and λ 5,6,7 can be complex.
In the CP-conserving limit, the neutral Higgs spectrum contains a CP-odd field A = S 3 and two CP-even scalars h and H which mix through the two-dimensional rotation matrix: We use the conventions M h ≤ M H , and 0 ≤α ≤ π so that sinα is always positive.

Yukawa sector
The 2HDM Yukawa sector is given by withΦ i (x) = iτ 2 Φ * i (x) the charge-conjugated scalar doublets with hypercharge Y = − 1 2 . Q L and L L denote the SM left-handed quark and lepton doublets, respectively, and u R , d R and R are the corresponding righthanded singlets, in the weak interaction basis.
The Yukawa couplings M f and Y f ( f = u, d, ) are complex 3 × 3 matrices which, in general, cannot be diagonalized simultaneously, generating FCNCs at tree level. This can be avoided by assuming that M f and Y f are proportional to each other [7]. In the masseigenstate fermion basis with diagonal matrices M f , one has then with arbitrary complex parameters ς f ( f = d, u, ), which introduce new sources of CP violation. The aligned 2HDM (A2HDM) Yukawa Lagrangian reads where P R,L ≡ 1±γ 5 2 , V is the CKM quark-mixing matrix and the neutral Yukawa couplings are given by The usual Z 2 symmetric models can be recovered with specific assignments of the alignment parameters.

Flavour misalignment
The alignment conditions (11) presumably hold at some high-energy scale Λ A and are spoiled by radiative corrections which induce a misalignment of the Yukawa matrices. However, the flavour symmetries of the A2HDM tightly constrain the possible FCNC structures, keeping their effects well below the present experimental bounds. The only FCNC local structures induced at one loop take the form [7,14], and absorbs the UV divergences from one-loop Higgspenguin diagrams in B 0 s,d → + − decays [1].

Effective Hamiltonian
The low-energy effective Hamiltonian describing B 0 s,d → + − decays is given by [15,16,17] where = e, µ, τ; q = d, s, and m b = m b (µ) denotes the b-quark MS running mass. Other possible operators are neglected because their contributions are either zero or proportional to the light-quark mass m q . The anomalous dimension of O 10 is zero due to the conservation of the (V − A) quark current in the massless quark limit. The operators O S and O P also have zero anomalous dimensions because the µ dependences of m b (µ) and the scalar current (qP R b)(µ) cancel each other. Therefore the Wilson coefficients C i do not receive additional renormalization from QCD corrections.

Calculation of the Wilson coefficients C 10,S,P
The Wilson coefficients C 10,S ,P are obtained by requiring the equality of one-particle irreducible amputated Green functions in the full and in the effective theories. The relevant Feynman diagrams for a given process can be created by the package FeynArts [18], with the model files provided by FeynRules [19]. The generated decay amplitudes are evaluated either with the help of FeynCalc [20], or using standard techniques such as the Feynman parametrization to combine propagators. We found full agreement between the results obtained with these two methods. Throughout the whole calculation, we set the light-quark masses m d,s to zero; while for m b , we keep it up to linear order.
In general, the Wilson coefficients C i are functions of the internal up-type quark masses, together with the corresponding CKM factors [21]: where x j = m 2 j /M 2 W , and F i (x j ) denote the loop functions. In deriving the effective Hamiltonian (17), the limit m u,c → 0 and the unitarity of the CKM matrix, have to be exploited. This implies that we need only to calculate explicitly the contributions from internal top quarks, while those from up and charm quarks are taken into account by means of simply omitting the massindependent terms in the basic functions F i (x t ).
The relevant Feynman diagrams are split into various box, penguin and self-energy diagrams, which are mediated by the top quark, gauge bosons, and Higgs scalars. In order to check the gauge independence of the final results, we perform the calculation both in the Feynman (ξ = 1) and in the unitary (ξ = ∞) gauges.

Wilson coefficients in the SM
In the SM, the dominant contribution to the decays B 0 s,d → + − comes from the Wilson coefficient C 10 , which arises from W-box and Z-penguin diagrams: where is the one-loop Inami-Lim function [22]. The factors η EW Y and η QCD Y account for the NLO electroweak [9] and the NNLO QCD corrections [10], respectively.
The coefficients C S and C P receive SM contributions from box, Z penguin, Goldstone-boson (GB) penguin and Higgs (h) penguin diagrams: The Goldstone contribution is of course absent in the unitary gauge. Explicit expressions can be found in [1].

Wilson coefficients in the A2HDM
In the A2HDM there are additional contributions from box and Z penguin diagrams, involving H ± exchanges, and from Higgs penguin diagrams. The only new contribution to C 10 comes from Z penguin diagrams and is gauge independent by itself: A2HDM 10 .
The Z penguin diagrams also generate contributions to C P . The sum of Z penguin diagrams and Goldstoneboson penguin diagrams is gauge independent: Neutral scalar exchanges induce both tree and loop diagrams. The loop contributions consist of the Higgspenguin and self-energy diagrams governed by the Yukawa couplings (12), whereas the tree ones are given by the misalignment couplings (15). The sum of these contributions can be written as: witĥ where λ The coefficients g 0 , g (a) j and g (b) j are functions of x t , x H + , ς u and ς d . g 0 and g (a) j are gauge independent because they do not involve any gauge bosons, while g (b) j are all related to Goldstone-boson vertices and, therefore, are identically zero in the unitary gauge. These g (b) j contributions cancel the gauge dependence from the box diagrams: The loop contributions with neutral scalar exchanges generate UV divergences which are cancelled by the renormalization of the misalignment coupling in (16). The µ dependence of the results is reabsorbed into the combination C R (M W ) = C R (µ) − ln (M W /µ).

Phenomenological analysis
Currently, only B s → µ + µ − is observed with a signal significance of ∼ 4.0 σ [13]. Thus we shall investigate the allowed parameter space of the A2HDM under the constraint from B(B 0 s → µ + µ − ). With updated input parameters, the SM prediction reads In order to explore constraints on the model parameters, it is convenient to introduce the ratio [15,16] where Γ s L(H) denote the lighter (heavier) eigenstate decay width of the B s meson, and ∆Γ s = Γ s L − Γ s H . The quantities S and P are defined as where the Wilson coefficients are given by: A2HDM 10 , Combining the SM prediction (32) with the latest experimental result (3), we get

Model parameters
We consider the CP-conserving limit and assume that the lightest CP-even scalar h corresponds to the observed neutral boson with M h 126 GeV. We have then 10 free parameters: 3 alignment couplings ς f , 3 scalar masses (M H , M A , M H ± ), 2 scalar-potential couplings (λ 3 , λ 7 ), the mixing angleα and the misalignment parameter C R (M W ). Four of them (α, λ 3,7 , C R (M W )) have minor impacts on R sµ , compared to the others. In order to simplify the analysis, we assign them the following values, using the bounds from earlier studies: The C S ,P contributions to (34) and (35) are suppressed by a factor M 2 B s /M 2 W . Therefore, unless there are large enhancements from the ς f parameters, the branching ratio shall be dominated by C 10 , where the A2HDM contribution depends only on |ς u | 2 and M H ± . We shall then discuss two possible scenarios: 1) |ς d, | |ς u | ≤ 2, where C 10 dominates, and 2) |ς d, | |ς u |, where C S and C P could play a significant role.

Small ς d,
The only relevant NP contribution is C A2HDM 10 which involves two parameters, ς u and M H ± . The constraints imposed byR sµ are shown in Fig. 1. In the left panel, we choose M H ± = 80, 200 and 500 GeV (upper, middle and lower curves, respectively). The shaded horizontal bands denote the allowed experimental region at 1σ (dark green), 2σ (green), and 3σ (light green), respectively. The right panel shows the resulting upper bounds on ς u , as function of M H ± . A 95% CL upper bound |ς u | ≤ 0.49 (0.97) is obtained for M H ± = 80 (500) GeV. Since C A2HDM 10 ∼ |ς u | 2 , this constraint is independent of any assumption about CP. For larger masses the constraint becomes weaker since the H ± contribution starts to decouple.

Large ς d,
In this case, C S and C P can induce a significant impact on the branching ratio. We vary ς d and ς within the range [−50, 50], and choose three representative values for ς u = 0, ±1. We also take three different representative sets of scalar masses: In Fig. 2, we show the allowed regions in the ς d -ς plane under the constraint fromR sµ . The regions with large ς d and ς are already excluded, especially when they have the same sign. The impact of ς u is significant: a nonzero ς u will exclude most of the regions allowed in the case with ς u = 0, and changing the sign of ς u Model − tan β cot β cot β Inert 0 0 0 will also flip that of ς . The allowed regions expand with increasing scalar masses, as expected, since the NP contributions gradually decouple from the SM.

2HDMs with discrete Z 2 symmetries
The usual Z 2 symmetric models are recovered for the values of ς f indicated in Table 1. In these models, the ratioR sµ only involves seven free parameters: M H ± , M H , M A , λ 3 , λ 7 , cosα and tan β. In the particular case of the type-II 2HDM at large tan β, our results agree with the ones calculated in Ref. [23]. It is also interesting to note that the B 0 s,d → + − branching ratios depend only on the charged-Higgs mass and tan β in this case. Fig. 3 shows the dependence ofR sµ on tan β, for three representative charged-Higgs masses: M H ± = 80, 200 and 500 GeV. The other two neutral scalar masses have been fixed at M H = M A = 500 GeV. The four different panels correspond to the models of types I, II, X and Y, respectively. A lower bound tan β > 1.6 is obtained at 95% CL under the constraint from the current experimental data onR sµ . This implies ς u = cot β < 0.63, which is stronger than the bounds obtained previously from other sources [14,24].

Conclusions
We have studied the rare decays B 0 s,d → + − within the general framework of the A2HDM. A complete oneloop calculation of the Wilson coefficients C 10 , C S and C P has been performed. They arise from various box, penguin and self-energy diagrams, as well as tree-level FCNC diagrams induced by the flavour misalignment interaction (15). The gauge independence of the results has been checked through separate calculations in the Feynman and unitary gauges, and the gauge relations among different diagrams have been examined in detail.
We have also investigated the impact of the current B(B 0 s → µ + µ − ) data on the model parameters, especially the resulting constraints on the three alignment couplings ς f . This information is complementary to the one obtained from collider physics and will be useful for future global data fits within the A2HDM.