$\bar{B}\to X_s \gamma$ in the Two Higgs Doublet Model up to Next-to-Next-to-Leading Order in QCD

We compute three-loop matching corrections to the Wilson coefficients $C_7$ and $C_8$ in the Two Higgs Doublet Model by applying expansions for small, intermediate and large charged Higgs boson masses. The results are used to evaluate the branching ratio of $\bar{B}\to X_s \gamma$ to next-to-next-to leading order accuracy, and to determine an updated lower limit on the charged Higgs boson mass. We find $\mhplus \ge 380 $GeV at 95% confidence level when the recently completed BABAR data analysis is taken into account. Our results for the charged Higgs contribution to the branching ratio exhibit considerably weaker sensitivity to the matching scale $\mu_0$, as compared to previous calculations.


Introduction
In view of missing (to date) New Physics signals at the Large Hadron Collider (LHC), it is of utmost importance to exploit precision calculations together with precise experimental results in order to look for deviations from the Standard Model (SM). In this context, the rare decayB → X s γ constitutes one of the most important processes. It is a loopgenerated Flavour-Changing-Neutral-Current (FCNC) transition, which makes it very sensitive to contributions from beyond-SM particles. Moreover, its branching ratio can be measured with an uncertainty of a few percent and, at the same time, the result can be predicted within perturbation theory with a similar uncertainty.
In this paper, we consider extensions of the SM Higgs sector by a second Higgs doublet, namely the so-called Two Higgs Doublet Models (2HDMs). They are constructed in such a way that no FCNC occur at the tree level [10]. Such models have five physical scalar degrees of freedom, among which there is a charged Higgs boson H ± that plays an important role forB → X s γ. We shall consider two versions of the model, usually denoted by 2HDM type-I and type-II where, respectively, either the same or two different Higgs doublet fields couple to the up-and down-type quarks. In these models, both Higgs doublets acquire vacuum expectation values v 1,2 such that v = v 2 1 + v 2 2 ≃ 246 GeV determines the W ± and Z boson masses in the same way as in the SM. The ratio v 2 /v 1 is denoted by tan β.
Comparison of the experimental results for B(B → X s γ) to predictions within the 2HDM type-II leads to the strongest constraint on the charged Higgs boson mass for tan β ∈ [1,25] (see, e.g., Sec. 5 of Ref. [11]; the precise range depends on the treatment of uncertainties). So far, the constraint has been derived using only the Next-to-Leading Order (NLO) expressions for the 2HDM contributions at the electroweak scale, while the SM contributions were treated at the NNLO [9]. In the present paper, we compute the missing two-and three-loop NNLO matching coefficients in the 2HDM, and re-do the analysis to extract a lower bound on M H + from the new experimental average (1).
The outline of the paper is as follows: In the next Section, we discuss the matching coefficients up to the three-loop order. Besides considering the 2HDM, we also re-compute the SM matching contribution, and improve the three-loop results for C 7 and C 8 . In Section 3, we use the new results to evaluate B(B → X s γ) to the NNLO accuracy, and to extract a lower bound on M H + . Section 4 contains our conclusions.

Formalism
The formalism to compute the Wilson coefficients in the 2HDM can be taken over from the SM case [12,13]. We shall follow the regularization and renormalization conventions of those papers. In particular, we adopt the effective Lagrangian and the definition of the operators P j (j = 1, . . . , 8,11) from Eqs. (2.1) and (2.2) of Ref. [13]. One should note that the dipole operators P 7 and P 8 are normalized there with inverse powers of the QCD coupling constant.
Since the additional degrees of freedom of the 2HDM are all heavy, they only influence the Wilson coefficients C i of the operators P i . In order to incorporate the 2HDM contribution in a manner that is analogous to the SM analysis of Ref. [13], we split the Wilson coefficients as where Q = c, t marks contributions from loop diagrams with virtual charm and top quarks, respectively.
Our matching calculation is performed at the renormalization scale µ 0 . It is chosen to be of the same order of magnitude as masses of the particles that are being decoupled (m t , M W and M H + ). For the NLO calculations within the SM we refer to Ref. [14]. The NNLO SM contributions to C Q,SM i (µ 0 ) have been computed in Refs. [12,13]. Suppression by m 2 c /M 2 W makes C c,2HDM i negligible. In the following, we shall consider C t,2HDM i only. It is convenient to decompose them as follows where C H(n) i is obtained from n-loop diagrams. For i = 7, 8, the tree-level coefficients C (µ 0 ) have been found in Refs. [15,16,17,18]. In the present paper, we compute the three-loop corrections C H(3) 7 (µ 0 ) and C H(3) 8 (µ 0 ), which constitutes the last missing element of the 2HDM Wilson coefficient evaluation that matters for B(B → X s γ) at the NNLO. Furthermore, we reproduce the two-loop corrections C H(2) i (µ 0 ) for i = 3, 4, 5, 6 that have been originally found in Ref. [19] and belong to the necessary NNLO matching in the 2HDM.
In the following two subsections, we discuss our results for the Wilson coefficients in the SM and 2HDM. Several variables turn out to be of convenience in this context:

Standard Model
Let us begin with recalling three-loop corrections to the matching coefficients in the SM. They depend on two masses: M W and m t . In Ref. [13], expansions of the three-loop diagrams for M W ≪ m t and M W ≈ m t have been performed. The final numerical results of that paper    [13]. In the present calculation, we have added eight more expansion terms for M W ≈ m t . They are shown in Fig. 1 as thin dashed lines and a solid curve that includes terms up to w 16 .
In the case of C t(3) 7 , we observe an overlap of the two expansions for 0.2 ≤ y ≤ 0.35, which gives us confidence that the exact curve is approximated with high accuracy by the Taylor expansion around w = 0 on one side, and by the asymptotic large-m t expansion on the other. The situation for C t(3) 8 in Fig. 1(b) is only slightly worse. We still observe an improvement w.r.t. Ref. [13] due to the additional terms in the w expansion. The vertical bands in Fig. 1 correspond to the current experimentally allowed region y = 0.492 ± 0.003 for µ 0 = m t . In the range 0.4 < y < 0.6, our improved results are very well approximated (to better than 0.1%) by which is consistent with Eq. (5) but much more accurate, and allows to substitute the updated value of y. In effect, the uncertainties get reduced by almost an order of magnitude.
As far as contributions from loops involving the charm quark are concerned, the corresponding coefficients in the range 0.4 < y < 0.6 can be obtained with high precision from the expressions given already in Ref. [13], namely They already provide a high-accuracy approximation in the physical region, so there is no need to consider higher-order terms in the expansions.

Two Higgs Doublet Model
In order to specify the notation, we provide the Lagrange density which defines interactions of the charged Higgs boson with fermions. Adopting the conventions from Ref. [15], we have where P L/R = (1 ∓ γ 5 )/2, V ij are the Cabibbo-Kobayashi-Maskawa matrix elements, u i and d j are the up-and down-type quarks with masses m u i and m d j , and G F is the Fermi constant. For the 2HDM of type-I and II, the couplings A d and A u take the values and respectively.
The 2HDM contributions to the Wilson coefficients in Eq. (2) are proportional to A i A ⋆ j . Since the terms involving A ⋆ d are suppressed by the strange-quark mass, we can safely neglect them. The remaining terms can be split into two parts as follows: For the computation of C t,2HDM i , we can use Eq. (5.1) of Ref. [13] that has been derived in the context of the SM. Its application to the 2HDM is straightforward after taking into account that, apart from the Wilson coefficients, also the electroweak counterterms (cf. Eqs. (4.7) to (4.10) of [13]) and the quantities B 7 and B 8 (cf. Eqs. (3.22) and (3.23) of [13]) receive additional contributions originating from the charged Higgs boson exchange. Decomposing each of these quantities as X = X SM + X 2HDM where X SM denotes the result given in Ref. [13], we obtain their 2HDM parts in D = 4 − 2ǫ dimensions by a simple one-loop calculation Both in the SM and 2HDM, the renormalization constants Z Q 0,sb and Z Q 2,sb enter the elec-troweak counterterm Lagrangian in the same way 1 All the remaining quantities appearing in Eq. (5.1) of Ref. [13] are precisely the same as in the SM.
In analogy to the SM, we have to consider vacuum integrals with two mass scales (M H + and m t ) in our matching calculation. Sample diagrams up to three-loops are shown in Fig. 2. At the one-and two-loop levels, the calculation can be performed exactly, and one obtains C 7 and C 8 as functions of m t /M H + [15,16,17,18]. At the three-loop level, we proceed as in Ref. [13], considering expansions around m t ≈ M H + , for m t ≪ M H + , and for m t ≫ M H + . In the first case, a simple Taylor expansion is sufficient, and we have computed the first 16 terms in u. Calculations involving strong hierarchies require non-trivial asymptotic expansions. In these cases, five terms in r and 1/r have been evaluated.
For the purpose of the present analysis, we have re-evaluated the Leading Order (LO) and NLO contributions to the renormalized Wilson coefficients, confirming the results of Refs. [15,16,17,18] and extending them to include higher powers in ǫ. Such an extension has been necessary for renormalization at the three-loop level. Two-loop (NNLO) Wilson coefficients of the four-quark-operators were calculated in Ref. [19] for the MSSM. We have performed an independent calculation of the charged Higgs boson contribution. Full agreement has been found. We refrain from displaying here explicit analytical results for the one-and two-loop Wilson coefficients, and refer to the Mathematica file available from Ref.
[20]. Let us only note that we set m t = m t (µ 0 ) in all those lower-order terms, which specifies our conventions for the NNLO expressions below.
At the three-loop level, we have been able to recover the dependence of C 7 and C 8 on µ 0 by applying renormalization group techniques to the analytical one-and two-loop results. We find In the corresponding Eq. (4.6) of Ref. [13], a factor of 4π e γǫ is missing in the global normalization, which we correct here. For clarity, the chirality projectors P R are now displayed explicitly.
Our results for C At this point, a comment concerning the expansions inū = 1 − m 2 t (µ 0 )/M 2 H + and u = 1 − M 2 H + /m 2 t (µ 0 ) is in order. For M H + ≈ m t , one can expand either inū or in u, and the two expansions are easily convertible into each other. 2 As it has been observed in Ref. [21], one expects better convergence of the expansion for M H + ≥ m t (M H + ≤ m t ) if the result is expressed in terms ofū (u). In the following, we shall always choose the better-suited representation without explicitly mentioning it.
In Figs. 3 and 4, we demonstrate that the above expansions are sufficient to obtain the Wilson coefficients for any M H + . In Fig. 3(a), the coefficient C   can define C H(3) 7,AuA ⋆ u (µ 0 = m t ) as a piecewise function using the expansions in the various limits. In Fig. 3(b) , the corresponding results for the coefficient C (µ 0 = m t ), which leads us to define

B(B → X s γ) in the 2HDM to NNLO
The framework for our numerical analysis is based on Ref. [8] where explicit results for the effective-theory description of B(B → X s γ) have been provided up to the NNLO. While the Wilson coefficients are known in a complete manner at this order, non-BLM NNLO corrections to the charm-quark-mass-dependent matrix elements (on-shell amplitudes) have been evaluated only in the large m c limit and extrapolated to the physical region.
Predictions within the 2HDM to be discussed below are obtained along the same algorithm and using the complete NNLO matching conditions from the previous Section. However, two-loop purely electroweak corrections to the matching are included in the SM part only, as they remain unknown in the 2HDM. One should keep in mind that such electroweak corrections and our new NNLO QCD matching ones may be of comparable size.
In the following, we shall discuss results for the 2HDMs of type-I and II that have been introduced in Eqs. (9) and (10). Most of the input parameters are adopted from Ref. [8], except for the strong coupling constant and the top quark mass for which we use the most up-to-date values that are given by [22,23,24] The corresponding MS top quark mass equals m t (m t ) = 163.5 GeV using three-loop accuracy in QCD [26,27] and neglecting electroweak effects. As in Ref. [8], our default value of In a first step, let us discuss the branching ratio dependence on tan β. In Fig. 5(a), we choose M H + = 400 GeV and show B(B → X s γ) in the 2HDM type-II for 0.5 ≤ tan β ≤ 10. The solid curve describes the NNLO result, while the dotted and dashed ones show the LO and NLO central values for comparison. One observes strong dependence for tan β < 2 and a nearly tan β-independent result for tan β > 2. Actually, from tan β = 10 to tan β = 50, the branching ratio changes only by 0.03%. In the following, tan β = 50 is going to be our default value for the type-II model; choosing tan β < 2 would strengthen the lower limit on M H + .
In Fig. 5 The main effect of our new three-loop terms is in reducing µ 0 -dependence of the decay rate and, in consequence, stabilizing the lower bound on M H + . This is illustrated in Fig. 6(a) where the charged Higgs contribution to the branching ratio ∆B ≡ B 2HDM −B SM is plotted as a function of µ 0 , while normalized to its own value at µ 0 = M H + = 400 GeV. Apart from the LO (dotted), NLO (dashed), NNLO (solid) curves, we also present the partial NNLO (dash-dotted) line that corresponds to the approach of Ref. [9]. Our calculation differs from the latter one precisely by including the 2HDM contributions to the NNLO matching. One observes a clear reduction of µ 0 -dependence when including higher order corrections. Whereas the partial NNLO result for the considered ratio varies by more than 6.6% when µ 0 is varied in the [80 GeV, 2M H + ] range, the corresponding variation of our present result with full NNLO matching remains below 1.6%. The overall size of ∆B(µ 0 ) amounts to around 25% of B(B → X s γ) SM in the considered case.
In order to determine a lower bound on M H + we follow the approach of Ref. [9], combining the experimental and theoretical uncertainties in quadrature. A one-sided 95% C.L. (99% C.L.) bound is obtained for the lowest value of M H + for which the difference between experimental and theoretical central values is 1.645 (2.326) times larger than the total uncertainty.
The theory uncertainty consists of four contributions that we take over from Ref. [8]. Let us briefly comment on each of them: • In Ref. [8], the non-perturbative uncertainty has been estimated to ±5%, which has been confirmed by the detailed investigation of Ref. [28]. We adopt this uncertainty for all the considered values of M H + .
• The charm quark mass dependence of the operator matrix elements is only partly known. Thus, an interpolation between the large-m c results [29] and reasonable assumptions for B(B → X s γ) at m c = 0 has to be performed. This uncertainty has been estimated in [8] to ±3%, which we again assume to be M H + -independent.
• The total parametric error is obtained by combining all the partial ones in quadrature. 4 It amounts to around 2 ÷ 3%, however, computed from scratch for each value of M H + .
• An estimate of higher order corrections is obtained by varying the renormalization scales in the ranges 80 GeV ≤ µ 0 ≤ max{2M H + , 320 GeV}, 1.25 GeV ≤ µ b ≤ 5 GeV and 1.224 This uncertainty is about 3 ÷ 4% but again we compute it for each value of M H + .
Contrary to Refs. [8,9], we work with asymmetric uncertainties resulting from the parameter and renormalization scale variation.
In Fig. 6(b), the NNLO result for B(B → X s γ) is shown as a function of M H + , together with an uncertainty band that is obtained by adding all the errors in quadrature. The dotted (M H + -independent) curves in Fig. 6(b)   From Fig. 6(b) one can extract (using the procedure described above) the following limits on M H + in the 2HDM type-II: M H + ≥ 380 GeV at 95% C.L. , M H + ≥ 289 GeV at 99% C.L. .
The above bounds replace the ones of Ref. [9] (295 and 230 GeV, respectively). A considerable improvement of the bounds arises mostly due to the shift in the experimental value of B(B → X s γ) that is now smaller than the one used in Ref. [9]. Our 95% C.L. limit is very close to the one presented in Ref. [7] (385 GeV) together with the new experimental average (Eq. (1)). On the other hand, it is significantly stronger than the one in Ref. [4] (327 GeV) that is based on the BABAR data alone. It is interesting to mention that when the matching scale µ 0 is varied between 80 and 400 GeV, the lower limits vary by around 25 GeV when our new three-loop 2HDM matching contributions are not included. This gets reduced to around 7 GeV only after including the three-loop corrections, which demonstrates the stabilizing effect of the full NNLO matching. On the other hand, for µ 0 = 160 GeV fixed, the correction strengthens the limit only slightly, by 5 ÷ 6 GeV.
The contour plots in Fig. 7   the projected uncertainty to be reached by Belle II and SuperB [30]. The current theory uncertainty has been used in Fig. 7(a), while Fig. 7(b) presents a future projection with assumed reduction of the theory uncertainty by a factor of two.
Let us now shortly discuss the 2HDM of type-I. In Fig. 8(a), we show the branching ratio for tan β = 1 as a function of M H + , including the uncertainties. The meaning of the curves is the same as in Fig. 6(b). Contrary to the type-II model, the branching ratio gets suppressed with respect to the SM, and the effect becomes larger for lower values of M H + , which implies increased discrepancy with respect to the experimental results. Thus, it is possible to set a lower bound on M H + . It is shown in Fig. 8(b) where the shaded area represents the part of the tan β-M H + plane that is excluded at 95% C.L. Similar results based on the partial NNLO predictions and the previous experimental average have been presented in Ref. [31].

Conclusions
Applying the method of expansion in mass ratios, we have evaluated three-loop matching conditions for the dipole operators P 7 and P 8 in the 2HDM. In effect, Wilson coefficients of all the operators that matter forB → X s γ in this model at the leading order in electroweak interactions are now known to the NNLO accuracy in QCD. This is true not only at the matching scale µ 0 but also at the low-energy scale µ b because the renormalization group evolution is the same as in the SM [32].
The main effect of including the NNLO matching is a significant reduction of µ 0 -dependence of the charged Higgs contribution to theB → X s γ branching ratio. It is particularly trans-parent when considering a lower bound on M H + in the type-II model. Before taking our correction into account, the bound varies by around 25 GeV when µ 0 is varied in a reasonable range [80, 400] GeV ∼ [ 1 2 m t , M H + ]. Including the correction reduces the variation by more than a factor of 3.
With the updated experimental average, we find that the 95% C.L. (99% C.L.) lower limit on M H + amounts to 380 (289) GeV in the 2HDM type-II. This is a universal (tan β-independent) bound that can only get stronger when tan β-dependence is taken into account. In practice, noticeable modifications occur for tan β smaller than around 2.
In the 2HDM type-I, a 95% C.L. lower limit on M H + from B(B → X s γ) can be derived for low tan β only, currently for tan β < 2.5 when M H + is above the LEP bound of around 80 GeV [23]. In this case, considerable reduction of µ 0 -dependence is observed, too.
With the semi-analytical results presented in this paper, constraints on the 2HDM at the NNLO level can easily be updated in the future, along with developments in the measurements and in calculations of the low-energy matrix elements (that are identical in the SM and in the 2HDM). More precise measurements are expected in a few years from Belle-II [33] and Super-B [34]. On the theoretical side, the main challenges are improvements in analyses of non-perturbative effects [28] together with perturbative calculations of the NNLO on-shell amplitudes beyond the large-m c limit. In the latter case, new results should become available soon [35].