Three-loop QCD corrections to B_s ->mu^+ mu^-

The decay B_s ->mu^+ mu^- in the Standard Model is generated by the well-known W-box and Z-penguin diagrams that give rise to an effective quark-lepton operator Q_A at low energies. We compute QCD corrections of order alpha_s^2 to its Wilson coefficient C_A. It requires performing three-loop matching between the full and effective theories. Including the new corrections makes C_A more stable with respect to the matching scale mu_0 at which the top-quark mass and alpha_s are renormalized. The corresponding uncertainty in |C_A|^2 gets reduced from around 1.8% to less than 0.2%. Our results are directly applicable to all the B_{s(d)} ->l^+ l^- decay modes.

Previous upper limits can be found in Refs. [5][6][7][8][9]. Although the experimental uncertainties are still quite large, they are expected to get significantly reduced within the next few years.
As far as the theory side is concerned, the B s meson decay into two muons is quite clean. In fact, the only relevant quantity that needs to be calculated at the leading order in α em and cannot be determined within perturbation theory is the leptonic decay constant f Bs . Its square enters the branching ratio as a multiplicative factor. Recent progress in the determination of f Bs from lattice calculations [10][11][12][13][14][15] gives a motivation for improving the perturbative ingredients, in particular the two-loop electroweak [16] and the three-loop QCD corrections.
Evaluation of the latter corrections is the main purpose of the present paper. Renormalization scale dependence of the truncated perturbation series is going to be significantly reduced. In our case, it refers to the branching ratio dependence on the scale µ 0 at which the top-quark mass and α s are renormalized. At the two-loop order, the corresponding uncertainty amounts to around 1.8%, which is a non-negligible component of the overall theoretical uncertainty.

(5)
In the SM, the operator Q A alone is sufficient because contributions from Q S and Q P to the branching ratio are suppressed by M 2 Bs /M 2 W with respect to that from Q A . In beyond-SM theories, the Wilson coefficients C S and C P can get enhanced, especially for an extended Higgs sector (see, e.g., Refs. [17,18]). Note that Q V = (bγ α γ 5 s)(μγ α µ) does not contribute at the leading order in α em due to the electromagnetic current conservation.
Using Eq. (3), the following result for the average time integrated branching ratio can be derived where Γ s H stands for the total width of the heavier mass eigenstate in the B sBs system. The quantities r, β and u are given by In the absence of beyond-SM sources of CP-violation, we have F P = 1 and F S = 1 − ∆Γ s /Γ s L , where Γ s L is the lighter eigenstate width, and ∆Γ s = Γ s L − Γ s H . In a generic case, from the results in Refs. [19,20] one derives where φ NP s describes the CP-violating "new physics" contribution to B sBs mixing, i.e. φ ccs s ≃ arg[(V * ts V tb ) 2 ] + φ NP s (see Sec. 2.2 of Ref. [21]).
In the SM, the branching ratio of B s → µ + µ − is proportional to the square of the Wilson coefficient C A which can be computed within perturbation theory. The calculation amounts to matching the amplitude 1 for s → b µ + µ − in the full SM to the one of the effective theory defined in Eq. (3). At the matching scale µ 0 , the W and Z bosons together with the top quark are integrated out simultaneously.
Barring higher-order electroweak (EW) corrections, the perturbative expansion of C A reads where α s ≡ α s (µ 0 ) in the MS scheme with five active quark flavours. No other definition of α s is going to be used throughout the paper. The one-loop term C A has been calculated 1 More precisely, we shall match thebsμµ one-light-particle-irreducible (1LPI) Green's functions at vanishing external momenta. for the first time in Ref. [22], and the two-loop correction C (1) A has been found in Refs. [23][24][25][26]. In this work, we compute the three-loop QCD correction C (2) A . Let us note that C (n) A are µ 0 -dependent, but C A itself is not, up to higher-order QED effects. It follows from the fact that the quark current in Q A is classically conserved in the limit of vanishing quark masses, while the chiral anomaly plays no role here, as we work at the leading order in flavour-changing interactions. Once the perturbation series on the r.h.s. of Eq. (9) is truncated, a residual µ 0 -dependence arises. Our present calculation aims at making this dependence practically negligible.
At each loop order, we shall split the coefficients C which are separately finite but gauge-dependent with respect to the EW gauge fixing.
Here, we use the background field version of the 't Hooft-Feynman gauge for the electroweak bosons, and the usual 't Hooft-Feynman gauge for the gluons. Most of the results have also been cross-checked using the general R ξ gauge for the gluons.
For the top quark mass renormalization, we shall always use the MS scheme in the full SM, i.e. m t ≡ m t (µ 0 ). The ratio m t /M W will enter our results via the following three variables: The ratio x is the only parameter on which the coefficients C (n) A depend, apart from the logarithms ln(µ 0 /M W ) or ln(µ 0 /m t ). For the Z-penguins, this is true after taking the leading-order EW relations between M Z , M W and sin 2 θ W into account.
Our paper is organized as follows: in the next two sections, we evaluate the matching coefficient C A up to three loops. Calculations of the W -boxes and the Z-penguins are discussed in Sections 2 and 3, respectively. Section 4 is devoted to a numerical analysis and examining the size of the evaluated three loop QCD corrections. We conclude in Section 5. Logarithmically enhanced QED corrections to C A are summarized in the Appendix.   Consequently, it is possible to write C W,(n) A in terms of the top-and charm-quark contributions where unitarity of the CKM matrix has been applied.
To obtain C W,t,(n) A and C W,c,(n) A , we compute off-shell 1LPI amplitudes both in the full theory and in the effective theory, and require that they agree at the scale µ 0 up to terms suppressed by heavy masses. In fact, on the full-theory side, all the external momenta can be set to zero, which leads to vacuum integrals up to three loops. On the other hand, in the effective theory, all loop corrections vanish in dimensional regularization after setting the external momenta and light quark masses to zero because the loop integrals are scaleless in this limit. Thus, we are only left with tree contributions.
There are basically two approaches to perform the matching. In the first one, the matching is performed in d = 4 − 2ǫ dimensions setting all the light masses strictly to zero. As a consequence, one generates spurious infrared divergences both in the full and effective theories. Such divergences cancel while extracting C A . However, due to the presence of additional poles in ǫ at intermediate steps, one has to introduce the so-called evanescent operators in the effective Lagrangian, which complicates the calculations. In an alternative matching procedure, finite light quark masses are introduced to obtain infrared and ultraviolet finite results, which allows for a matching in four dimensions. In the latter case, no evanescent operators matter. In the following, we describe both matching procedures in more detail.

Matching in d dimensions
The evanescent operator which enters the effective Lagrangian when the matching is performed in d dimensions reads [25]  Note that this operator vanishes in d = 4 dimensions, and thus the limit d → 4 can only be taken after matching in d dimensions.
Before performing the matching, we have to replace the combination C A Q A + C E A Q E A by the corresponding renormalized expression that can be written as [25] where Z ψ is the MS quark wave function renormalization constant. Loop corrections to Z ψ , Z N N , Z EE and Z N E contain no finite parts 2 but at most poles in ǫ. As far as Z EN is concerned, we require that amplitudes proportional to C E A vanish for d → 4. In consequence, Z EN may contain both pole parts and (uniquely defined) finite terms. For our purpose, the renormalization constants are needed up to two loops.
The renormalization constants Z N N , Z N E , Z EN and Z EE are computed from the diagrams like those in Fig. 2, with insertions of Q A and Q E A . Since we are only interested in ultraviolet poles of momentum integrals, all the masses can be set to zero, and an external momentum q flowing through the quark lines is introduced. Our results read where n f = 5 denotes the number of active quark flavours. The results for Z N N and Z N E are true to all orders in QCD due to the (already mentioned) quark current conservation for massless quarks. Concerning Z EN , we confirm the one-loop result from Ref. [25], whereas the two-loop expression is new. Note that Z EE does not matter for our calculation, and thus we have left its two-loop part unevaluated.
In the first step of our matching calculation, we determine the s → bµ + µ − transition amplitude in the full theory, where the Dirac structure of each Feynman diagram is projected onto Q A and Q E A (see, e.g., Ref. [27]). This gives us the unrenormalized amplitudes, which we denote by C W A,bare and C E A,bare , respectively. In this step, vacuum diagrams up to three loops with two different mass scales have to be computed. Although some classes of Feynman diagrams of this type have been studied in the literature (see, e.g., Ref. [28]), we have decided to perform expansions in various limits, which leads to handy results for the matching coefficients. Actually, we follow the same strategy as in Refs. [29,30], namely, we expand in the limits M W ≪ m t and M W ≈ m t , i.e. y ≪ 1 and w ≪ 1, where terms up to order y 12 and w 16 are evaluated, respectively. A simple combination of the two expansions provides an approximation to the three-loop contribution, which for all practical purposes is equivalent to an exact result.
The actual calculation has been performed with the help of QGRAF [31] to generate the Feynman diagrams, q2e and exp [32] for the asymptotic expansions [33] and MATAD [34], written in Form [35], for evaluation of the three-loop diagrams. We have performed our calculation for an arbitrary gauge parameter in QCD, and have checked that it drops out in our final result for the matching coefficient.
For renormalization of the full-theory contributions, we need the one-loop renormalization constant for the QCD gauge coupling Here, N ǫ = (µ 2 0 /m 2 t ) ǫ e γǫ Γ(1 + ǫ) makes the renormalized α s in the full SM equal to the MS-renormalized α s in the five-flavour effective theory, to all orders in ǫ. As far as the top quark mass is concerned, its two-loop MS renormalization constant Z mt in the full SM is expressed in terms of the above-defined α s , which gives Furthermore, for the wave-function renormalization, only the difference between the renormalization constants in the full and effective theories has to be taken into account (see Section 4 of Ref. [29]): At this point, all the ingredients are available to perform the matching according to the following equations (Q = c, t): where ∆T E,t, (1) and ∆T W,t,(n) denote contributions from the top-quark mass counterterms which can be written as Here, the following notation has been used: "m bare t → Z mt m t " means that the bare top quark mass is replaced by the renormalized one times the renormalization constant. Afterwards, we expand in α s and take the coefficient at [α s /(4π)] n (n = 1, 2), which is indicated by the subscript at the round bracket.
Our final results for the evanescent Wilson coefficients up to two loops read The results for C W,t A and C W,c A will be given in Subsection 2.4.

Matching in four dimensions
In order to have a cross check of the results for C W A from the previous subsection, we have performed the matching also for infrared finite quantities, which can be done in four dimensions avoiding evanescent operators [25]. No spurious infrared divergences arise when small but non-vanishing masses are introduced for the strange and bottom quarks.
In the full theory, this leads to Feynman diagrams with up to four different mass scales. We evaluate them using asymptotic expansions in the limit In addition, we use either m t ≫ M W or m t ≈ M W , as in the previous subsection. All the external momenta are still set to zero. Asymptotic expansions are conveniently performed with the help of exp [32].
On the effective-theory side, the loop corrections do not vanish any more due to the finite quark masses. We compute the necessary one-and two-loop Feynman integrals in the limit After renormalization of the two-loop expression on the effective-theory side and the three-loop result on the full-theory side, the finite parts are matched for ǫ → 0. After the matching, it is possible to take the limit m s → 0 and m b → 0. This way, we obtain the same results as in the previous calculation where the infrared divergences have been regulated using dimensional regularization.
Although we only had to compute the leading non-vanishing contributions in the light quark masses, the calculational effort has been significantly higher than for the matching in d dimensions described in the previous subsection. Thus, we have applied the method with light masses only to cross check the first two (three) terms in the expansion in y (w), using a general R ξ gauge though.

Results
At the one-and two-loop orders, we have confirmed the results with full dependence on x from Ref. [25], and evaluated in addition terms up to O(ǫ 2 ) and O(ǫ), respectively. For completeness, we present the results for ǫ → 0 which are given by Analytic expressions including O(ǫ) terms can be downloaded from [36].
With the help of the exact two-loop result, we can extract the full x-dependence in front of the ln µ 0 terms at the three-loop level. We find The corresponding coefficient in the top sector is given by are accurate to better than 1% in the corresponding quantities.

General remarks
The second type of contribution to C A arises from the so-called Z-boson penguins. Sample diagrams at the one-, two-and three-loop orders are shown in Fig. 4. In contrast to the W -box diagrams, there is no contribution from evanescent operators to C Z,(n) A . However, flavour non-diagonal loop contributions to the light quark kinetic terms require introduction of an EW counterterm which already appears at the one-loop level. The corresponding counterterm Lagrangian (see Eq. (13) of Ref. [30] 3 ) can be written in the following form where D µ is the covariant derivative involving the neutral gauge boson fields (Z, γ, g). While only the one-loop contributions to Z Q 2,sb were needed in theB → X s γ case [29], now also the two-and three-loop corrections of order α s and α 2 s matter. The two-loop ones were also necessary in Refs. [23][24][25][26][27]. Perturbative expansions of Z Q 2,sb are conveniently written as For determination of Z Q 2,sb , a two-point function with incoming strange quark and outgoing bottom quark has to be considered. Sample diagrams at one, two and three loops are shown in Fig. 5. We refrain from explicitly listing the results but refer to [36] for computerreadable expressions.  The counterterm Z Q 2,sb is either inserted in the tree-level amplitude or in two-loop diagrams containing a closed top quark loop on the gluon propagator, as shown in Fig. 6. Insertions of the counterterm into other loop contributions lead to massless tadpoles which vanish in dimensional regularization.

Fermion triangle contribution
There is a class of Feynman diagrams which require special attention, namely those containing a closed triangle quark loop (see Fig. 7). For these contributions, a naive treatment of γ 5 as anticommuting is not possible, and a more careful investigation is necessary. We have followed two approaches which are described below. Similarly to the anomaly cancellation in the SM, contributions with the up, down, strange and charm quarks running in the triangle loop cancel pairwise within each family. Thus, only the top and bottom quarks need to be considered, as the top is the only massive quark in our calculation.
In our first approach, we adopt the prescription from Ref. [37] and replace the axial-vector coupling in the triangle loop as follows In a next step, we pull out the ε tensor and take the trace of the loop diagram in d dimensions. In the resulting object, we perform the replacements and proceed from now on in the same way as with the other diagrams contributing to C Z A . In the second approach, we do not take the trace in the triangle loop at all, but only use the cyclicity property for traces and anticommutation relations for the γ matrices (not for γ 5 ) in order to put γ 5 to the end of each product under the trace. Afterwards, we perform the tensor loop integration, and use again anticommutation relations to bring the resulting expressions to the form where only the axial-vector part of the (Z boson)-lepton coupling has been taken into account. In a next step, we add and subtract 24γ µ ⊗ γ µ γ 5 to the first, and 24γ µ γ 5 ⊗ γ µ γ 5 to the second structure in Eq. (35). This way, we obtain the Wilson coefficients for the trace evanescent operators [38] and a contribution to C A . Actually, the latter is given by (−24) times the prefactor of the second structure in Eq. (35).
The two methods, which lead to identical results for C Z A , have been applied both to the three-loop diagrams themselves and to the counterterm contributions (cf. Fig. 7).

Matching formula
In analogy to Eq. (12) we can write where C Z,Q,tria.
A are the contributions from the triangle diagrams described in the previous subsection.
The calculation of C Z,(n) A proceeds along the same lines as for the W -box contribution. In particular, we set all the external momenta to zero, and expand the Feynman integrals in the full theory both for m t ≫ M W and m t ≈ M W . Furthermore, we renormalize the topquark mass, α s and the wave function in analogy to the W -box case. As before, all loop corrections vanish in the effective theory, which finally leads to the following matching K Q,(n) denote tree-level contributions from the EW counterterm (31) which take a simple form Finally,K Q stands for the counterterm contributions from two-loop diagrams like those in Fig. 6.
We observe that after inserting explicit results on the right-hand side of Eq. (38), all the terms proportional to sin 2 θ W cancel out, and C Z,t A becomes independent of the weak mixing angle. This can be understood by recalling similarities between the Z boson and the photon couplings to other particles in the background field gauge, as well as the structure of the counterterm in Eq. (31). Once the quark kinetic terms in the effective theory are imposed to be flavour diagonal, the same must be true for dimension-four quarkphoton couplings. In effect, the counterterm in Eq. (31) automatically renormalizes away all the zero-momentum quark-(Z boson) interactions that come with sin 2 θ W .
Another interesting thing to note is that C Z,c,(n) A = 0 at each loop order in the background field version of the 't Hooft-Feynman gauge, which means that only the triangle contributions are non-vanishing in the charm sector. One of the ways to understand this fact is again by considering diagrams where the Z boson (together with the muons) is replaced by an external photon.

Results
With our calculation, we have confirmed the one-and two-loop results from Ref. [23] which are given by (for ǫ → 0) Furthermore, similarly to C W A , we obtain exact dependence on M W and m t for the µ 0dependent terms which read For the generic three-loop contributions, terms up to order y 12 and w 16 have been evaluated, as in the W -box case in Section 2. In a numerical form, they read In the limit of large top quark mass, the coefficient C Z,(2) A grows as m 2 t , which has its origin in the Yukawa interaction of the charged pseudo-goldstones with the top quark. For this reason, we plot in Fig. 8 the combination y 2 C Z,(2) A where sums of the results from Eqs. (43) and (44) are shown as dashed (M W ≪ m t ) and solid lines (M W ≈ m t ). Note that after multiplication by y 2 , the latter is expanded in w = 1 − y 2 . Again, one observes that the two approximations coincide for y ≈ 0.4, which suggests that a combination of the two expansions covers the whole range between y = 0 and y = 1. In the physical region, the expansion around M W = m t provides an excellent approximation. It is interesting to note that the fermion triangle contribution is more than an order of magnitude larger than C Z,t,(2) A . This is particularly true for the physical value of y where we have y 2 C Z,t,(2) A ≈ −0.02.
A handy approximation formula which works to better than 1% for 0.3 < y < 0.7 reads

Numerical analysis
In this section, we shall discuss numerical effects of our three-loop QCD corrections. The B(B s → µ + µ − ) branching ratio in the SM is proportional to |C A | 2 (cf. Eq. (6) with C S = C P = 0 and F P = 1). Here, we shall consider |C A | 2 only. Evaluation of the branching ratio itself is relegated to a parallel article [39] where also the new two-loop EW corrections [16] are included.
The relevant parameters are as follows. For the gauge boson masses, we take M Z = 91.1876 GeV [40] and M W = 80.358 GeV (calculated from G F , M Z and α em ). For the strong coupling, α s (M Z ) = 0.1184 in the five-flavour QCD is used [40]. Four-loop renormalization group equations (RGE) are applied to evolve it to other scales. For the topquark mass, our input is M t = 173.1 GeV [40] which we treat as the pole mass. We convert it to the MS scheme with respect to QCD, but include no shift due to the EW interactions. This means that our m t should be understood as renormalized on-shell with respect to the EW interactions. As far as QCD is concerned, we use a three-loop relation for converting M t to m t (m t ), which gives m t (m t ) ≃ 163.5 GeV. Next, four-loop RGE are used to find m t (µ 0 ) at other values of µ 0 . Fig. 9 shows the matching scale dependence of |C A | 2 . The dotted, dashed and solid curves show the leading order (LO), next-to-leading-order (NLO) and next-to-next-to-leadingorder (NNLO) results, respectively. In the current case, they correspond to one-, twoand three-loop matching calculations.
The range of the plot corresponds roughly to µ 0 ∈ 1 2 M W , 2m t , which might be considered reasonable given that both the W -boson and the top quark are decoupled simultaneously. However, our Wilson coefficient has a trivial RGE (at the LO in EW interactions) but it is quite sensitive to m t . In consequence, the main reason for its µ 0 -dependence here is the top-quark mass renormalization. Thus, for estimating uncertainties due to truncation of the QCD perturbation series at each order, we shall use a more narrow range µ 0 ∈ 1 2 m t , 2m t . One observes in Fig. 9 that the prediction for |C A | 2 has already improved a lot after including the NLO QCD corrections. The variation at the NLO level amounts to around 1.8% only, for µ 0 ∈ 1 2 m t , 2m t . Once the new three-loop corrections are taken into account, the uncertainty gets reduced to less than 0.2%, which can be treated as negligible for all practical purposes.
Another conclusion that can be drawn from Fig. 9 is that for µ 0 = m t the NLO correction is moderate (2.2%), while the NNLO correction essentially vanishes. Although µ 0 = m t has been anticipated to be an optimal scale in the past [24], there has been no convincing theoretical argument for such a choice. Our explicit three-loop calculation has been actually necessary to suppress the QCD matching uncertainties in |C A | 2 to the current sub-percent level.
For µ 0 = 160 GeV, our final result for C A is well approximated by the fit which is accurate to better than 0.1% for α s (M Z ) ∈ [0.11, 0.13] and M t ∈ [170, 175] GeV.
Let us stress that the O(α em ) term in Eq. (46) stands for both the NLO EW matching corrections at µ = µ 0 , as well as effects of the evolution of C A down to µ = µ b ∼ m b , according to the RGE. Once the QCD logarithms get resummed, the latter effects behave not only like O(α em ), but also like O(α em /α s ) and O(α em /α 2 s ), which means that they are potentially more important than the NNLO QCD corrections evaluated here. However, the actual numerical impact of the O(α em /α s ) and O(α em /α 2 s ) terms on the decay rate amounts to around −1.5% only [41], which has been checked using anomalous dimension matrices from Refs. [42,43]. The necessary expressions are given in the Appendix.
As far as the NLO EW matching corrections are concerned, they have been known for a long time only in the m t ≫ M W limit [44]. A complete calculation of these two-loop corrections has recently been finalized [16]. Their numerical effect depends on the scheme used at the LO. A detailed discussion of this issue is presented in Ref. [16]. Let us only mention that the semi-perfect stabilization of µ 0 -dependence in Fig. 9 at the NNLO in QCD takes place only because we have renormalized m t and M W on shell with respect to the EW interactions. If we used MS at µ 0 for the EW renormalization of m t and M W , then acceptable stability would be observed only after including the very two-loop EW corrections.

Conclusions
We have evaluated the NNLO QCD corrections to the Wilson coefficient C A that parametrizes the B s → µ + µ − branching ratio in the SM. For this purpose, three-loop matching between the SM and the relevant effective theory has been performed. Tadpole integrals depending on m t and M W have been evaluated with the help of expansions starting from the limits m t ≈ M W and m t ≫ M W , which for all practical purposes is equivalent to an exact calculation. When masses of the light quarks and their momenta are set to zero, care has to be taken about the so-called evanescent operators, similarly to the NLO case [25]. Such operators have also been helpful in dealing with diagrams where γ 5 was present under traces.
Our results for the renormalized matching coefficients C W A and C Z A can be downloaded in a computer-readable form from [36]. Including the new corrections makes C A more stable with respect to the matching scale µ 0 at which the top-quark mass and α s are renormalized. Apart from B s → µ + µ − , our calculation is directly applicable to all the B s(d) → ℓ + ℓ − decay modes, and it matters for other processes mediated by Z-penguins and W -boxes, e.g.,B → X s νν, K → πνν, or short-distance contributions to K L → µ + µ − . However, it is only B s → µ + µ − for which the three-loop accuracy is relevant at present.
An updated SM prediction for B(B s → µ + µ − ) is presented in a parallel article [39] where also the new two-loop EW corrections [16] are included. and the Graduiertenkolleg "Elementarteilchenphysik bei höchster Energie und höchster Präzision". M.M. acknowledges partial support by the National Science Centre (Poland) research project, decision DEC-2011/01/B/ST2/00438.

Appendix: Logarithmically enhanced QED corrections
In this appendix, we present explicit expressions for the logarithmically-enhanced QED corrections to C A . Beyond the LO in α em , its perturbative expansion at the matching scale µ 0 reads where C s A stands for the scale-independent O(α 0 em ) contribution as given in Eq. (9). Using the RGE from Refs. [42,43] one obtains the following result at the scale µ b where G includes all the NLO EW corrections that are not logarithmically enhanced. The quantities F 1,2,3 depend on η = αs(µ 0 ) αs(µ b ) and x = m 2 t M 2 W . We find with z ≃ 0.0553 and the remaining coefficients summarized in Table 1. The functions Y (x), V (x) and E(x) originate from the one-loop SM matching conditions [22] for various operators in the effective theory. They read Y (x) = 3x 2 8(x − 1) 2 ln x + The coefficients in Table 1 satisfy the following identities: With the help of them one can easily check that the terms in Eq. (48) proportional to F 1,2,3 are finite in the limit α s → 0. Logarithms ln n µ 2 0 µ 2 b with n = 1, 2 arise in this limit, which explains why the corresponding QED corrections are called logarithmically enhanced. The QED logarithms α em ln in the corresponding terms that brings inverse powers of α s into the final results.