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

The rare decays $B_{s,d}^0 \to \ell^+\ell^-$ are analyzed within the general framework of the aligned two-Higgs doublet model. We present a complete one-loop calculation of the relevant short-distance Wilson coefficients, giving a detailed technical summary of our results and comparing them with previous calculations performed in particular limits or approximations. We investigate the impact of various model parameters on the branching ratios and study the phenomenological constraints imposed by present data.


Introduction
The recent discovery of a Higgs-like boson [1,2], with properties compatible with the Standard Model (SM) expectations [3][4][5], is one of the greatest achievements in the past decades in particle physics and represents a major confirmation of our present theoretical paradigm. The LHC data suggest that the electroweak symmetry breaking (EWSB) is probably realized in the most elegant and simple way, i.e., via the Higgs mechanism implemented through one scalar SU(2) L doublet. An obvious question we are now facing is whether the discovered 126 GeV state corresponds to the unique Higgs boson incorporated in the SM, or it is just the first signal of a much richer scenario of EWSB. None of the fundamental principles of the SM forbids the possibility of an enlarged scalar sector associated with the EWSB.
Among the many possible scenarios for new physics (NP) beyond the SM, the two-Higgs doublet model (2HDM) [6] provides a minimal extension of the scalar sector that naturally accommodates the electroweak (EW) precision tests, giving rise at the same time to a large variety of interesting phenomenological effects [7]. The scalar spectrum of the model consists of two charged fields, H ± , and three neutral ones, h, H and A, one of which is to be identified with the Higgs-like boson found at the LHC. The direct search for these additional scalar states at high-energy collisions, or through indirect constraints via precision flavour experiments, is an important task for the next years. This will also be helpful to gain further insights into the scalar sector of supersymmetry (SUSY) and other models with similar scalar contents.
Within the SM, flavour-changing neutral current (FCNC) interactions are forbidden at tree level, and highly suppressed at higher orders, due to the Glashow-Iliopoulos-Maiani (GIM) mechanism [8]. In a generic 2HDM, however, tree-level FCNC interactions generally exist, through non-diagonal couplings of neutral scalars to fermions. The unwanted FCNCs can be eliminated, imposing on the Lagrangian an ad-hoc discrete Z 2 symmetry; depending on the different possible Z 2 charge assignments, this results in four types of 2HDMs (I, II, X and Y) [7], all satisfying the hypothesis of natural flavour conservation (NFC) [9]. A more general alternative is to assume the alignment in flavour space of the Yukawa matrices for each type of right-handed fermions [10]. The so-called aligned two-Higgs doublet model (A2HDM) results in a very specific structure, with all fermion-scalar interactions being proportional to the corresponding fermion masses. It also contains as particular cases the different versions of the 2HDM with NFC, while at the same time introduces new sources of CP violation beyond the Cabibbo-Kobayashi-Maskawa (CKM) phase [11]. These features make the A2HDM a very interesting theoretical framework, which leads to a rich and viable phenomenology, both in high-energy collider experiments [12,13], as well as in low-energy flavour physics [14,15].
In the field of rare B-meson decays, the purely leptonic processes B 0 s,d → + − , with = e, µ or τ , play an outstanding role in testing the SM and probing physics beyond it, because they are very sensitive to the mechanism of quark-flavour mixing. Within the SM, the FCNC transition is mediated by a one-loop amplitude, suffers from a helicity-suppression factor m /m b , and is characterized by a purely leptonic final state. The first two features result in a double suppression mechanism, responsible for the extremely rare nature of these decays. The third feature implies that these processes are theoretically very clean, with the only hadronic uncertainty coming from the B-meson decay constants f B s,d . All these considerations make the rare leptonic decays B 0 s,d → + − a formidable probe of physics beyond the SM, especially of models with a non-standard Higgs sector like multi-Higgs doublet models [16][17][18][19][20] as well as various SUSY scenarios [17,18,[20][21][22].
As far as the experimental side is concerned, the decay modes with = µ are especially interesting because the corresponding final state can be easily tagged. Over the last decade the upper bounds for the branching ratios of these decays have been improving continuously, thanks to the CDF and DØ collaborations at the Tevatron and, more recently, the ATLAS, CMS and LHCb experiments at the LHC [23]. In November 2012, the LHCb experiment reported the first evidence of the decay B 0 s → µ + µ − , at the 3.5 σ level [24]. The signal significance has been raised, respectively, to 4.0 σ and 4.3 σ by LHCb and CMS, after analyzing the currently available data set, with the averaged time-integrated branching ratio given by 3.0 +1.0 −0.9 × 10 −9 CMS [26] , where the CMS uncertainty includes both the statistical and systematic components, but is dominated by the statistical uncertainties. The two measurements lead to the weighted world average [27] B(B 0 s → µ + µ − ) exp. = (2.9 ± 0.7) × 10 −9 .
At the same time, the branching fraction of B 0 d → µ + µ − has also been determined with a signal significance of 2 σ by the two experiments: 3.5 +2.1 −1.8 × 10 −10 CMS [26] .
The corresponding combined result reads [27] These measurements are in remarkable agreement with the latest updated predictions within the SM [28]: where the next-to-leading order (NLO) corrections of EW origin [29], as well as the QCD corrections up to the next-to-next-to-leading order (NNLO) [30], have been taken into account.
Although the experimental uncertainties are still quite large, they are expected to get significantly reduced within the next few years [31]. All these experimental and theoretical progresses will lead to new stringent constraints on physics beyond the SM.
Motivated by the above considerations, in this work we shall perform a study of the rare leptonic decays B 0 s,d → + − within the A2HDM. Our paper is organized as follows. In section 2 we give a brief overview of the A2HDM Lagrangian, especially of its Yukawa and scalar sectors.
In section 3 we summarize the SM results and describe the full one-loop calculation of the relevant Feynman diagrams in the A2HDM. We have performed the calculation in two different gauges, Feynman (ξ = 1) and unitary (ξ = ∞), in order to check the gauge-independence of our results. In section 4 we discuss the impact of the model parameters on the branching ratios of these decays, taking into account the latest implications from the LHC Higgs data. Our conclusions are made in section 5. Finally, the appendix contains the explicit results for the individual Higgs-penguin diagrams.

The aligned two-Higgs doublet model
The 2HDM extends the SM with the addition of a second scalar doublet of hypercharge Y = 1 2 [6]. In the so-called "Higgs basis", in which only one doublet gets a nonzero vacuum expectation value, the two doublets can be parametrized as where G ± and G 0 denote the Goldstone fields, and v = ( √ 2G F ) −1/2 246 GeV. The five physical scalar degrees of freedom are given by the two charged fields H ± (x) and three neu- The latter are related with the S i fields through an orthogonal transformation, which is fixed by the scalar potential: are 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 right-handed singlets, in the weak interaction basis. All fermionic fields are written as 3-vectors in flavour space and, accordingly, the couplings M f and Y f (f = u, d, ) are 3 × 3 complex matrices.
In general, the Yukawa matrices M f and Y f cannot be simultaneously diagonalized in flavour space. Thus, in the fermion mass-eigenstate basis with diagonal mass matrices M f , the corresponding Yukawa matrices Y f remain non-diagonal, giving rise to tree-level FCNC interactions.
In the A2HDM, the tree-level FCNCs are eliminated by requiring the alignment in flavour space of the two Yukawa matrices coupling to a given type of right-handed fermions [10] where the three proportionality parameters ς f (f = d, u, ) are arbitrary complex numbers and introduce new sources of CP violation. The Yukawa interactions of the physical scalars with the fermion mass-eigenstate fields then read [10] where P R,L ≡ 1±γ 5 2 are the right-handed and left-handed chirality projectors, M f the diagonal fermion mass matrices, and V the CKM quark-mixing matrix [11]. The couplings of the neutral scalar fields to fermion pairs are given by In the A2HDM, all fermionic couplings to scalars are proportional to the corresponding fermion masses, and the only source of flavour-changing interactions is the CKM quark-mixing matrix V , while all leptonic couplings and the quark neutral-current interactions are diagonal in flavour. All possible freedom allowed by the alignment conditions is encoded by the three family-universal complex parameters ς f , which provide new sources of CP violation without tree-level FCNCs [10]. The usual models with NFC, based on discrete Z 2 symmetries, are recovered for particular values of the couplings ς f , as indicated in Table 1. Explicit examples of symmetry-protected underlying theories leading to a low-energy A2HDM structure have been discussed in Ref. [33].
The alignment conditions in Eq. (17) presumably hold at some high-energy scale Λ A and are spoiled by radiative corrections. These higher-order contributions induce a misalignment of the Yukawa matrices, generating small FCNC effects suppressed by the corresponding loop factors [10,14,34,35]. However, the flavour symmetries of the A2HDM tightly constrain the possible FCNC structures, keeping their effects well below the present experimental bounds [14,15]. Using the renormalization-group equations (RGEs) [35], one can check that the only FCNC local structures induced at one loop take the form [14,34] which vanishes identically when ς d = ς u (Z 2 models of types I, X and inert) or ς d = −1/ς * u (types II and Y), as it should be.
Although the numerical effect of the local term in Eq. (20) is suppressed by m q m 2 q /v 3 and quark-mixing factors, its tree-level contribution is needed to render finite the contribution from one-loop Higgs-penguin diagrams to B 0 s,d → + − , as will be detailed later. The renormalization of the coupling constant C is determined to be where D is the space-time dimension. Thus, the renormalized coupling satisfies Assuming the alignment to be exact at the scale Λ A , i.e., C R (Λ A ) = 0, this implies C R (µ) = ln (Λ A /µ).

Effective Hamiltonian
The rare leptonic B 0 s,d → + − decays proceed through loop diagrams in both the SM and the A2HDM. After decoupling the heavy degrees of freedom, including the top quark, the weak gauge bosons, as well as the charged and neutral Higgs bosons, these decays are described by a low-energy effective Hamiltonian [36][37][38] where G F is the Fermi coupling constant, α = e 2 /4π the QED fine-structure constant, and s W = sin θ W the sine of the weak angle. The effective four-fermion operators are given, respectively, In the SM, the contributions from the scalar and pseudoscalar operators are quite suppressed and, therefore, are usually neglected in phenomenological analyses. However, they can be much more sizeable in models with enlarged Higgs sectors, such as the A2HDM, especially when the Yukawa and/or scalar-potential couplings are large. Therefore, the B 0 s,d → + − data provide useful constraints on the model parameters. To get the theoretical predictions for B(B 0 s,d → + − ), the main task is then to calculate the three Wilson coefficients C 10,S,P in both the SM and the A2HDM, details of which will be presented in the next few subsections.

Computational method
The standard way to find the Wilson coefficients is to require equality of one-particle irreducible amputated Green functions calculated in the full and in the effective theory [39]. The former requires the calculation of various box, penguin and self-energy diagrams. We firstly use the program FeynArts [40], with the model files provided by the package FeynRules [41], to generate all the Feynman diagrams contributing to the decays B 0 s,d → + − , as well as the corresponding amplitudes, which can then be evaluated straightforwardly.
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. As the external momenta are much smaller than the masses of internal top-quark, gauge bosons, as well as charged and neutral scalars, the Feynman integrands are expanded in external momenta before performing the loop integration [42] where k denotes the loop momentum, M a heavy mass and l an arbitrary external momentum.
In addition, we employ the naive dimensional regularization scheme with an anti-commuting γ 5 to regularize the divergences appearing in Feynman integrals. After the Taylor expansion and factorizing out the external momenta, the integrals remain dependent only on the loop momentum and the heavy masses M . Subsequently, we apply the partial fraction decomposition [43] 1 which allows a reduction of all the Feynman integrals to those in which only a single mass parameter occurs in the propagator denominators. Finally, after reduction of tensor integrals to scalar ones, the only non-vanishing one-loop integrals take the form [44] with an arbitrary integer power n and with m = 0.
The computational procedure has also been checked through an independent analytic calculation of the Feynman diagrams, using more standard techniques such as the Feynman parametrization to combine propagators. We found full agreement between the results obtained with these two methods.
It should be noted that, in deriving the effective Hamiltonian in Eq. (23), the limit m u,c → 0 and the unitarity of the CKM matrix, have been implicitly exploited. In general, the Wilson coefficients C i are functions of the internal up-type quark masses, together with the corresponding CKM factors [39]: where x j = m 2 j /M 2 W , and F i (x j ) denote the loop functions. The unitarity relation in Eq. (28) implies vanishing coefficients C i if the internal quark masses are set to be equal, i.e., x u = x c = x t . For this reason, 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 mass-independent terms in the basic functions F i (x t ). For simplicity, we also introduce the following mass ratios: where m t = m t (µ) is the top-quark running mass in the MS scheme, and h SM the SM Higgs boson.
In order to make a detailed presentation of our results, we shall split the different contributions to the Wilson coefficients into the form: The pieces labeled with "SM" only involve SM fields (without the Higgs), while those denoted In B 0 s,d → + − the external momenta are small compared to the EW scale M W . One can then set all external momenta to zero when evaluating C 10 . However, the external momenta must be taken into account to evaluate the scalar Wilson coefficients C S and C P , otherwise some contributions would be missed.

Wilson coefficients in the SM
In the SM, the dominant contributions to the decays B 0 s,d → + − come from the W -box and Z-penguin diagrams shown in Figs. 1 and 2, respectively, which generate the Wilson coefficient: where  is the one-loop function that was calculated for the first time in Ref. [47]. The factor η EW Y accounts for both the NLO EW matching corrections [29], as well as the logarithmically enhanced QED corrections that originate from the renormalization group evolution [28,30], while the coefficient η QCD Y stands for the NLO [48,49] and NNLO [30] QCD corrections.
When the small external momenta are taken into account, the SM W -box and Z-penguin diagrams also generate contributions to the Wilson coefficients C S and C P . The contribution from diagram 1.2 can be neglected, because it contains two leptonic Goldstone couplings which generate a suppression factor m 2 /M 2 W . The scalar contribution from the remaining box diagrams is given by: where the two different expressions correspond to the results obtained in the Feynman and unitary gauges, respectively.
In the SM there is an additional contribution to the scalar Wilson coefficient C S from the Higgs-penguin diagrams shown in Fig. 3, which is by itself gauge dependent [46,50,51] and should cancel the gauge dependence of the W -box contribution. We find the result: The sum of the two contributions to C S is indeed gauge independent: The contribution from the SM W -box diagrams ( Fig. 1) to the pseudoscalar Wilson coefficient C P is given by: Additional contributions to C P are generated by the Z-and Goldstone-penguin diagrams shown in Figs. 2 and 4, respectively. The contributions from diagrams 4.6, 4.7 and 4.8 are proportional to the light-quark mass and can be therefore neglected. We find: C GB penguin, SM and Using the above results, one can easily check that the SM contribution to C P is also gauge independent: The GIM mechanism has eliminated those contributions which are independent of the virtual top-quark mass. However, the ln x t terms in the Wilson coefficients C SM S and C SM P do not vanish in the massless limit: at These infrared-sensitive terms arise from diagrams 1.1 and 2.1 in both gauges. The corresponding contributions from virtual up and charm quarks cancel in the matching process with the lowenergy effective theory, which has the same infrared behaviour. 2

Wilson coefficients in the A2HDM
In the A2HDM, the only new contribution to C 10 comes from the Z-penguin diagrams shown in Fig. 5. The result is gauge independent and given by In the particular case of the type-II 2HDM (or MSSM), ς u = 1/ tan β, this result agrees with the one calculated in Ref. [20].
The box diagrams shown in Fig. 6 involve charged scalar exchanges and contribute to the Wilson coefficients C A2HDM S and C A2HDM P . The contributions from diagrams 6.3 and 6.4 can be 2 In the low-energy effective theory the same ln x c (ln x u ) terms appear from analogous diagrams with a cν (uν ) or cc (uū) loop connecting two four-fermion operators. neglected, since they are proportional to m 2 /M 2 W . For the scalar coefficients we find the results: while the pseudoscalar contributions are given by: Most of the previous calculations in the literature focused on the type-II 2HDM in the large tan β limit; i.e., only those contributions proportional to tan 2 β were kept, which correspond to the ς d ς * terms in Eqs. (48)- (51). For this specific case, our results agree with Ref. [16].
Similarly to the SM case, the coefficient C A2HDM P receives additional contributions from Zand Goldstone-penguin diagrams shown in Figs. 5 and 7, respectively. They are given by: C GB penguin, A2HDM The gauge dependence of these two contributions compensates each other. Since there is no contribution from Goldstone-penguin topologies in the unitary gauge, the Z-penguin result should satisfy in this case: This relation has been validated by the actual calculation.

Neutral scalar exchange
The Wilson coefficients C A2HDM S and C A2HDM P receive a direct tree-level contribution from the scalar-exchange diagram shown in Fig. 8, where the FCNC vertex ϕ 0 is b is generated by the local operator in Eq. (20). This contribution must be combined together with the scalar penguin diagrams shown in Fig. 9. The structure of the common ϕ 0 i¯ vertex relates the resulting scalar and pseudoscalar Wilson coefficients, which take the form: The contributions from diagrams 9.4, 9.7, 9.8 and 9.14 are proportional to the light-quark mass m q and, therefore, vanish in our massless approximation. Diagrams 9.1, 9.3, 9.11 and 9.13 in Feynman gauge and diagrams 9.1, 9.3, 9.5, 9.6, 9.9 and 9.10 in unitary gauge generate a divergent contribution, which is not eliminated by the GIM mechanism; i.e., it remains even after summing over contributions of the three virtual up-type quarks. This divergence matches exactly the expected behaviour predicted through the RGEs, which originated in the local term L FCNC . Thus, the one-loop divergence is cancelled by the renormalization of the coupling C in Eq. (21) which, moreover, reabsorbs the µ dependence of the loops into the combination The scalar penguin diagrams 9.2, 9.12, 9.15 and 9.16 involve the cubic couplings ϕ 0 i H + H − , ϕ 0 i G + G − , ϕ 0 i H + G − and ϕ 0 i G + H − , respectively, which are functions of the scalar-potential parameters. Since the last three couplings can be fully determined in terms of the vacuum expectation value v and the scalar masses and mixings, we can express the total scalar-exchange (tree-level plus one-loop penguin) contribution in the form: where λ j (x t , x H + , ς u , ς d ) and g limit ς u,d → 0, x H,A → ∞, x h → x h SM , R i2,i3 → 0, R 11 → 1, this result reduces to the SM expression in Eqs. (38) and (39).

The orthogonality relation [13]
allows us to separate the total contribution from the functions g j (x t , x H + , ς u , ς d ), which does not depend on the neutral scalar masses: It is also noted that the functions g j (x t , x H + , ς u , ς d ) only receive contributions in the Feynman gauge, because they arise from the scalar penguin diagrams involving the Goldstone bosons.
Actually, the gauge dependent pieces from the box diagrams shown in Figs. 1 and 6 are exactly cancelled by these terms: The remaining contributions in Eq. (56), which are all proportional to 1/M 2 ϕ 0 i , are gauge independent but are sensitive to the scalar mixing parameters. Nevertheless, a naive mixingindependent estimate can be obtained in the limit of degenerate neutral-scalar masses: We shall perform our phenomenological analyses in the CP-conserving limit, with real potential and alignment parameters, where A = S 3 is a CP-odd state while H and h are two CP-even states defined by the rotation in Eq. (11). The 1/x ϕ 0 i contributions take then the form: where cα = cosα and sα = sinα. For degenerate neutral scalars, this reproduces the results in Eqs. (63) and (64) (in the CP-conserving limit).
The terms proportional to C R (M W ) in Eqs. (65) and (66) are absent in Z 2 -symmetric models, because the alignment conditions are protected by the Z 2 symmetry at any scale. In the particular case of the type-II 2HDM at large tan β, the only terms enhanced by a factor tan 2 β originate from the ς g (for C S ) and ς g 3 (see Eqs. (112) and (113)). In this specific case, our results agree with the ones calculated in Ref. [16]. Especially, we confirmed the observation that the dependence on the masses of the neutral Higgs bosons from the penguin and fermion self-energy diagrams drops out in their sum without invoking any relation between the mixing angle and the Higgs masses [16].

B 0 s,d → + − branching ratio
Due to the pseudoscalar nature of the B q meson, only the following two hadronic matrix elements are involved in B 0 s,d → + − decays: where f Bq and M Bq are the B q -meson decay constant and mass, respectively. The second equation follows from the first one by using the QCD equation of motion for the quark fields.
Starting with Eq. (23) and using Eq. (67), we can express the branching ratio of B 0 s,d → + − decays as where τ Bq is the B q -meson mean lifetime, and P and S are defined, respectively, as [36,37] P ≡ We have approximated the negligibly small (and usually neglected) SM scalar/pseudoscalar In the SM, P = 1 and S = 0. In a generic case, however, P and S can carry nontrivial CP-violating phases φ P and φ S . It is also noted that, even in models with comparable Wilson coefficients, the contributions from O S and O P are suppressed by a factor M 2 Bq /M 2 W with respect to that from O 10 . Therefore, unless there were large enhancements for C S and C P , the coefficient C 10 still provides the dominant contribution to the branching ratio.
In order to compare with the experimental measurement, the effect of B 0 q −B 0 q oscillations should be taken into account, and the resulting averaged time-integrated branching ratio is given by [36,37] where A ∆Γ is a time-dependent observable introduced firstly in Ref. [37], and y q is related to the decay width difference ∆Γ q between the two B q -meson mass eigenstates, with Γ q H(L) denoting the heavier (lighter) eigenstate decay width and Γ q = τ −1 Bq the average B q -meson width. Within the SM, A ∆Γ = 1 and the averaged time-integrated branching ratio is given by (73) 3 Here, C SM S and C SM P denote the full SM contribution, including the Higgs-penguin terms.
where the second line is valid only in the absence of beyond-SM sources of CP violation, which will be assumed in the following. 4 4 Numerical results

Input parameters
To evaluate numerically the branching ratios in Eqs. (73) and (74), we need several input parameters collected in Table 2. For the matching scale µ 0 ∼ O(M W ) and the low-energy scale µ b ∼ O(m b ), we fix them to µ 0 = 160 GeV and µ b = 5 GeV [28]. In addition, the onshell scheme is adopted for the EW parameters, which means that the Z-boson and top-quark masses coincide with their pole masses, and the weak angle is given by where M W = 80.359 ± 0.012 GeV is the W -boson on-shell mass obtained according to the fit formulae in Eqs. (6) and (9) of Ref. [59].
For the top-quark mass, we assume that the combined measurement of Tevatron and LHC [53] corresponds to the pole mass, but increase its systematic error by 1 GeV to account for the intrinsic ambiguity in the m t definition; i.e. we shall take M t = (173.34 ± 0.27 ± 1.71) GeV.
With the aid of the Mathematica package RunDec The decay constants f Bq are taken from the updated FLAG [54] average of N f = 2+1 lattice determinations, which are consistent with the naive weighted average of N f = 2 + 1 [61][62][63] and N f = 2 + 1 + 1 [64,65] results. For the B q -meson lifetimes, while a sizable decay width difference ∆Γ s has been established [55], the approximation 1/Γ d H 1/Γ d L ≡ τ B d can be safely set, given the tiny SM expectation for ∆Γ d /Γ d [66].
For the CKM matrix element |V cb |, we adopt the recent inclusive fit performed by taking into account both the semileptonic data and the precise quark mass determinations from flavourconserving processes [56]. However, one should be aware of the present disagreement between inclusive and exclusive determinations [54]. With |V cb | fixed in this way, the needed CKM factors are then obtained (within the SM) from the accurately known ratio |V * tb V ts /V cb | [57, 58].

SM predictions
Within the SM, only the Wilson coefficient C SM 10 is relevant and, using the fitting formula in Eq. (4) of Ref. [28] (which has been transformed to our convention for the effective Hamiltonian), The EW and QCD factors introduced in Eq. (34) are extracted as: With the input parameters collected in Table 2 where a 1.5% nonparametric uncertainty has been set to the branching ratios, and the main parametric uncertainties come from f Bq and the CKM matrix elements [28]. The small differences with respect to the results given in Ref. [28] are due to our slightly different (more conservative) input value for the top-quark mass M t .
In order to explore constraints on the model parameters, it is convenient to introduce the ratio [36,37] where the hadronic factors and CKM matrix elements cancel out. Combining the theoretical SM predictions in Eq. (77) with the experimental results in Eqs. (2) and (4), we get to be compared with the SM expectation R SM sµ = R SM dµ = 1. Since only the B s → µ + µ − branching ratio is currently measured with a signal significance of ∼ 4.0σ [27], we shall investigate the allowed parameter space of the A2HDM under the constraint from R sµ given in Eq. (79). Although the experimental uncertainty is still quite large, it has already started to put stringent constraints on many models beyond the SM [36].
Notice that, in addition to modifying the ratios R q , the scalar contributions to B 0 q -B 0 q mixings also change the fitted values of the relevant CKM parameters and, therefore, the normalization B(B 0 q → + − ) SM . This should be taken into account, once more precise B 0 q → + − data becomes available, through a combined global fit.

Choice of model parameters
In the following we assume that the Lagrangian of the scalar sector preserves the CP symmetry i.e., that the only source of CP violation is still due to the CKM matrix. This makes all the alignment and scalar-potential parameters real. Assuming further that the lightest CPeven scalar h corresponds to the observed neutral boson with M h 126 GeV, there are ten free parameters in our calculation: three alignment parameters ς f , three scalar masses (M H , , one mixing angleα, two scalar-potential couplings (λ 3 , λ 7 ), and the misalignment parameter C R (M W ).
In order to gain insight into the parameter space allowed by B 0 s,d → + − decays, it is necessary to take into account information about the h(126) collider data and flavour physics constraints, as well as EW precision observables, which will be crucial for making simplifying assumptions and reducing the number of relevant variables. Explicitly, the following constraints and assumptions on the model parameters are taken into account: • Firstly, the mixing angleα is constrained at | cosα| > 0.90 (68% CL) through a global fit to the latest LHC and Tevatron data for the h(126) boson [13], which is very close to the SM limit; i.e., the lightest CP-even scalar h behaves like the SM Higgs boson.
• To assure the validity of perturbative unitarity in the scalar-scalar scattering amplitudes, upper bounds on the quartic Higgs self-couplings are usually imposed by requiring them to be smaller than 8π [7]; i.e., |λ 3,7 | 8π. • The charged Higgs mass is assumed to lie in the range M H ± ∈ [80, 500] GeV, which would require |ς u | ≤ 2 to be compatible with the present data on loop-induced processes, such as Z →bb, b → sγ and B 0 s,d −B 0 s,d mixing, as well as the h(126) decays [13,15].
• The alignment parameters ς d and ς are only mildly constrained through phenomenological requirements involving other model parameters. As in our previous works, we restrict them at |ς d, | ≤ 50 [15].
• At present, there are no useful constraints on the misalignment parameter C R (M W ). For simplicity, it is assumed to be zero.
Numerically, it is found that the ratio R sµ is less sensitive to the scalar-potential couplings λ 3 and λ 7 than to the other model parameters, especially when the alignment parameters are small and/or the neutral scalar masses are large. The mixing angleα, when constrained in the range cosα ∈ [0.9, 1], is also found to have only a marginal impact on R sµ . Thus, for simplicity, we shall assign the following values to these parameters: As can be seen from Eqs. (69) and (70), the Wilson coefficients C S and C P are always accompanied with the power-suppressed factor M 2 Bq /M 2 W compared to C 10 . The NP contribution to C 10 is, however, proportional to |ς u | 2 and depends only on the charged-scalar mass. It is, therefore, interesting to discuss the following two special cases with respect to the choice of the alignment parameters: The first one is when |ς d, | |ς u | ≤ 2, where the NP contribution is dominated by C 10 while C S and C P are negligible. The second one is when |ς d, | |ς u |, which means that C S and C P play a significant role.

Small ς d,
When the alignment parameters ς d, are of the same size as (or smaller than) ς u , the NP contributions from C S and C P are negligible. In this case, we need only to focus on the Wilson coefficient C 10 , which is the sum of the SM contribution C SM 10 and the charged-Higgs contribution C A2HDM 10 due to Z-penguin diagrams shown in Fig. 5. The latter involves only two free parameters, ς u and M H ± , and goes to zero when ς u → 0 and/or M H ± → ∞.
The dependence ofR sµ on the alignment parameter ς u with three typical charged-Higgs masses (80, 200 and 500 GeV) is shown in Fig. 10. One can see that, with the contributions from C S and C P ignored, the observableR sµ puts a strong constraint on the parameter ς u .
For M H ± = 80 (500) GeV, a 95% CL upper bound |ς u | ≤ 0.49 (0.97) is obtained, with the assumption |ς d, | |ς u |, which is stronger than the constraint from R b [14]. Since C A2HDM 10 ∼ |ς u | 2 , this constraint is independent of any assumption about CP and, therefore, applies in the most general case. 5 For larger charged-Higgs masses, the constraint becomes weaker as the NP effect starts to decouple, reflected by lim which cover the lower, intermediate, and upper range, respectively, of the allowed scalar spectrum.
With the above specification, we show in Fig. 11 the allowed regions in the ς d -ς plane under the constraint fromR sµ . One can see that, irrespective of the scalar masses, regions with large ς d and ς are already excluded, especially when they have the same sign. The impact of ς u , even when varied within the small range [−1, 1], is found to be significant: a nonzero ς u will exclude most of the regions allowed in the case with ς u = 0, and changing the sign of ς u will also flip that of ς . This is mainly due to the factors ς 2 d ς * u appearing in the functions g (a) 2 and g (a) 3 defined, respectively, by Eqs. (112) and (113). It is also observed that the allowed regions expand with increasing scalar masses, as expected, since larger scalar masses make the NP contributions gradually decouple from the SM.

Z 2 symmetric models
The five types of Z 2 -symmetric models listed in Table 1 are particular cases of the CP-conserving A2HDM, with the three alignment factors ς f reduced to a single parameter tan β = v 2 /v 1 ≥ 0.
In the particular scalar basis where the discrete Z 2 symmetry is implemented, the scalarpotential couplings µ i and λ i must be real, and µ 3 = λ 6 = λ 7 = 0; however, the rotation into the Higgs basis generates non-zero values of µ 3 = − 1 2 λ 6 v 2 and λ 7 . Furthermore, the alignment condition is protected by the Z 2 symmetry at any energy scale, which means that the misalignment parameter C R (M W ) does not contribute and the Higgs-penguin diagrams are free of divergences. Thus, for Z 2 -symmetric models, the ratioR sµ only involves seven free parameters: M H ± , M H , M A , λ 3 , λ 7 , cosα, and tan β.
A much more constrained case is the inert 2HDM, where the Z 2 symmetry is imposed in the Higgs basis: all SM fields and Φ 1 are even while Φ 2 → −Φ 2 under the Z 2 transformation.
This implies that there is no mixing between the CP-even neutral states h and H, and the scalars H, A and H ± decouple from the fermions: cosα = 1, λ 6 = λ 7 = 0, ς f = 0. Moreover, the couplings of h to fermions and vector bosons are identical to the SM ones. Therefore, in the inert modelR inert sµ = 1.
For the other four types of Z 2 -symmetric models, we continue to use the assignments cosα = 0.95 and λ 3 = λ 7 = 1. One can easily check that the effects of M H and M A onR sµ are tiny, unless tan β is extremely small which is excluded by the flavour constraint |ς u | ≤ 2. For simplicity, we fix them to be M H = M A = 500 GeV in the following analysis.  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,15].
One can see that the predictedR sµ in the type-I, type-X and type-Y models are almost indistinguishable from each other and, in the large tan β region, approach the SM prediction, irrespective of the choices of scalar masses. For the type-II model, on the other hand, an enhancement ofR sµ is still possible in the large tan β region. This can be understood since the Wilson coefficients in the type-II model contain the factor tan 2 β arising from the product of alignment parameters ς f , while in the other three models they contain at most one power of tan β. So only the type-II model can receive a large tan β enhancement, which has been studied intensively in the literature [16][17][18]. It is also interesting to note that in the type-II 2HDM with large tan β the B 0 s,d → + − branching ratios depend only on the charged-Higgs mass and tan β [16].

Conclusions
In this paper, we have performed a detailed analysis of the rare decays B 0 s,d → + − within the general framework of the A2HDM. Firstly, we presented a complete one-loop calculation of the short-distance Wilson coefficients C 10 , C S and C P , which arise from various box and penguin diagrams, and made a detailed technical summary of our results and a comparison with previous calculations performed in particular limits or approximations. In order to make sure our results are gauge independent, the calculations were carried out in both the Feynman and the unitary gauges.
With the current data on B(B 0 s → µ + µ − ) taken into account, we have also investigated the impact of various model parameters on the branching ratios and studied the phenomenological constraints imposed by present data. The resulting information about the model parameters will be crucial for the model building and is complementary to the collider searches for new scalar resonances in the near future.
When |ς d, | |ς u |, the contributions to B(B 0 s → µ + µ − ) from the scalar and pseudoscalar operators are negligible compared to the leading Wilson coefficient C 10 . Since C A2HDM 10 ∼ |ς u | 2 , the measured B(B 0 s → µ + µ − ) branching ratio implies then an upper bound on the up-family alignment parameter, which only depends on the charged Higgs mass. At 95% CL, we obtain: for M H ± = 80 (500) GeV and |ς d, | |ς u | .
This bound is stronger than the constraints obtained previously from other sources [14,15].
The role of the scalar and pseudoscalar operators becomes much more important for large values of |ς d, |. This region of parameter space was previously explored within the context of the type-II 2HDM, where these contributions are enhanced by a factor tan 2 β. Our analysis agrees with previous results in the type-II case and shows, moreover, that this tan 2 β enhancement is absent in the Z 2 -symmetric models of types I, X and Y, which approach the SM prediction for large values of tan β. From the current experimental data onR sµ , we derive the 95% CL bound: tan β > 1.6 , for 2HDMs of types I, II, X and Y.
This implies ς u = cot β < 0.63, which is also stronger than the bounds obtained previously from other sources [14,15]. It would be interesting to analyze the possible impact of the new CP-violating phases present within the A2HDM framework, at large values of |ς d, |. They could generate sizeable phases φ P and φ S in Eqs. (69) and (70), which could manifest themselves in the time-dependence of the B 0 s → µ + µ − decay amplitude [36]. To quantify the possible size of this effect requires a more careful assessment of the allowed parameter space of the A2HDM, which we plan to further investigate in future works.
In the unitary gauge, we find: In the Feynman gauge the results are: