NNLO QCD counterterm contributions to $B$ \to $X_s \gamma$ for the physical value of $m_c$

TL;DR

This paper calculates the ultraviolet counterterm contributions at the physical charm quark mass to improve the precision of the NNLO QCD corrections in the $B o X_s \gamma$ decay rate, reducing theoretical uncertainties.

Contribution

It provides the first explicit calculation of all relevant ultraviolet counterterms for the physical charm mass, refining previous interpolation-based estimates.

Findings

Enhanced accuracy in the NNLO QCD correction estimates.

Reduced theoretical uncertainties in the $B o X_s \gamma$ branching ratio.

Improved understanding of charm quark mass effects in rare B decays.

Abstract

One of the most important O($\alpha_s^2$) corrections to the $B \to X_s \gamma$ branching ratio originates from interference of contributions from the current-current and photonic dipole operators. Its value has been estimated using an interpolation in the charm quark mass between the known results at $m_c=0$ and for $m_c \ll m_b/2$. An explicit calculation for the physical value of $m_c$ is necessary to remove the associated uncertainty. In the present work, we evaluate all the ultraviolet counterterm contributions that are relevant for this purpose.

Tables

Table 1. Table 1: Master integrals with unitarity cuts indicated by the dashed lines. Thick solid, thin solid and dotted internal lines denote propagators with masses m b subscript 𝑚 𝑏 m_{b} , m c subscript 𝑚 𝑐 m_{c} and zero, respectively. Thick solid external lines come with the momentum p 𝑝 p such that p 2 = m b 2 superscript 𝑝 2 superscript subscript 𝑚 𝑏 2 p^{2}=m_{b}^{2} . Dots on the propagators indicate that they are squared. Arrows in I 16 subscript 𝐼 16 I_{16} define an irreducible numerator (see Eq. ( 3.2 )).

$I_1$				$I_7$				$I_{13}$
$I_2$				$I_8$				$I_{14}$
$I_3$				$I_9$				$I_{15}$
$I_4$				$I_{10}$				$I_{16}$
$I_5$				$I_{11}$				$I_{17}$
$I_6$				$I_{12}$				$I_{18}$

Equations

$${\mathcal{B}}_{s\gamma }^{\mathrm{SM}}=(3.36\pm 0.23)\times 10^{-4},$$

$${\mathcal{B}}_{s\gamma }^{\exp }=(3.32\pm 0.15)\times 10^{-4}.$$

$${\mathcal{B}}_{s\gamma }^{\exp }=(3.27\pm 0.14)\times 10^{-4}$$

$$Q_1=({\bar{s}}_L{\gamma }_{\mu }{T}^a{c}_L)({\bar{c}}_L{\gamma }^{\mu }{T}^a{b}_L),Q_2=({\bar{s}}_L{\gamma }_{\mu }{c}_L)({\bar{c}}_L{\gamma }^{\mu }{b}_L),Q_7=\frac{e{m}_b}{16{\pi }^2}({\bar{s}}_L{\sigma }^{\mu \nu }{b}_R)F_{\mu \nu }.$$

$$\Gamma (b\to X_s^p\gamma {)}_{E_{\gamma }>E_0}=\frac{G_F^2{m}_b^5{\alpha }_{\mathrm{e}m}}{32{\pi }^4}{\left|V_{ts}^{\ast }{V}_{tb}\right|}^2\sum _{i,j=1}^8{C}_i({\mu }_b)C_j({\mu }_b){\hat{G}}_{ij},$$

$${\hat{G}}_{ij}={\hat{G}}_{ij}^{(0)}+{\tilde{\alpha }}_{\mathrm{s}}{\hat{G}}_{ij}^{(1)}+{\tilde{\alpha }}_{\mathrm{s}}^2{\hat{G}}_{ij}^{(2)}+\mathcal{O}({\tilde{\alpha }}_{\mathrm{s}}^3).$$

$${\tilde{\alpha }}_{\mathrm{s}}{\hat{G}}_{i7}^{(1)}+{\tilde{\alpha }}_{\mathrm{s}}^2{\hat{G}}_{i7}^{(2)}$$

$${\hat{G}}_{i7}^{(2)\mathrm{bare}},{\hat{G}}_{i7}^{(1)\mathrm{bare}},{\hat{G}}_{i7}^{(1)3P},{\hat{G}}_{i7}^{(1)m},\text{mbox}fori=1,2,$$

$${\hat{G}}_{7j}^{(1)\mathrm{bare}},\text{mbox}forj=11,12.$$

$${\hat{G}}_{17}^{(1)\mathrm{bare}}=-\frac{1}{6}{\hat{G}}_{27}^{(1)\mathrm{bare}},{\hat{G}}_{17}^{(1)3P}=-\frac{1}{6}{\hat{G}}_{27}^{(1)3P},{\hat{G}}_{17}^{(1)m}=-\frac{1}{6}{\hat{G}}_{27}^{(1)m},{\hat{G}}_{7(11)}^{(1)\mathrm{bare}}=-\frac{1}{6}{\hat{G}}_{7(12)}^{(1)\mathrm{bare}}.$$

$$Q_{11}$$

$$Q_{12}$$

$${\hat{G}}_{7(12)}^{(1)\mathrm{bare}}=-4ϵ\left\{(1+ϵ){\hat{G}}_{27}^{(1)2P(d)}+(5+ϵ)\left[{\hat{G}}_{27}^{(1)2P(u)}+{\hat{G}}_{27}^{(1)3P}\right]\right\},$$

$${\hat{G}}_{27}^{(1)\mathrm{bare}}$$

$${\hat{G}}_{27}^{(1)m}$$

$${\hat{G}}_{27}^{(1)2P(d)}$$

$${\hat{G}}_{27}^{(1)m,2P(d)}$$

$$I_{13}$$

$$I_{16}$$

$${\mu }^{2ϵ}\int d{\mathrm{PS}}_2{W}_2$$

$${\mu }^{4ϵ}\int d{\mathrm{PS}}_3{W}_3$$

$${\hat{G}}_{77}^{(0)}=8\pi (1-ϵ){\mu }^{2ϵ}\int d{\mathrm{PS}}_2=\frac{e^{\gamma ϵ}\Gamma (2-ϵ)}{\Gamma (2-2ϵ)}.$$

$$I_1^{\prime }=\frac{2-2ϵ}{z}{I}_1,I_2^{\prime }=\frac{1-ϵ}{z}{I}_2,I_3^{\prime }=\frac{1-ϵ}{z}{I}_3,I_4^{\prime }=\frac{2-2ϵ}{z(1-4z)}(I_1-I_2)-\frac{2-4ϵ}{1-4z}{I}_4,$$

$$I_5^{\prime }=-2I_6,I_6^{\prime }=\frac{(1-ϵ{)}^2}{z^2(1-4z)}{I}_1-\frac{2-7ϵ+6{ϵ}^2}{z(1-4z)}{I}_5-\frac{(6-16ϵ)z+ϵ}{z(1-4z)}{I}_6,$$

$$I_7^{\prime }=-2I_8,I_8^{\prime }=\frac{(1-ϵ{)}^2}{4z^2(1-z)}(I_1-2I_2)-\frac{(1-2ϵ)(2-3ϵ)}{4z(1-z)}{I}_7+\frac{1-4ϵ-(3-8ϵ)z}{2z(1-z)}{I}_8,$$

$$I_9^{\prime }=\frac{1-ϵ}{z}(I_9-I_4)+\frac{2-3ϵ}{z}{I}_7+2I_8,I_{10}^{\prime }=-\frac{1}{z}{I}_6+\frac{1-2ϵ}{z}{I}_{10},$$

$$I_{11}^{\prime }=\frac{(1-ϵ{)}^2}{(1-2ϵ)z^2}{I}_1+\frac{1-ϵ}{z}{I}_4+\frac{1-3ϵ}{z}{I}_{11},I_{12}^{\prime }=\frac{(1-ϵ{)}^2}{2(1-2ϵ)z^2}{I}_1+\frac{1-ϵ}{z}{I}_4-\frac{2-3ϵ}{2z}{I}_7-I_8-\frac{ϵ}{z}{I}_{12}-I_{13},$$

$$I_{13}^{\prime }=\frac{(1-ϵ{)}^2}{2z^2(1-z)(1-4z)}(3I_2-\frac{(2-3z+4z^2)}{2z}{I}_1)+\frac{2-6ϵ+4{ϵ}^2}{z(1-4z)}{I}_4-\frac{(1-2ϵ)(2-3ϵ)}{4z(1-z)}{I}_7-\frac{(1-2ϵ)(1+z)}{2z(1-z)}{I}_8-\frac{2ϵ}{z}{I}_{13},$$

$$I_{14}^{\prime }=\frac{1-ϵ}{z}{I}_{14},I_{15}^{\prime }=\frac{-2+2ϵ}{z(1-4z)}{I}_{14}-\frac{2-2ϵ+(2-4ϵ)z}{z(1-4z)}{I}_{15}-\frac{6-8ϵ}{z(1-4z)}{I}_{16},I_{16}^{\prime }=\frac{1-ϵ}{z}{I}_{15}+\frac{3-4ϵ}{z}{I}_{16},$$

$$I_{17}^{\prime }=-\frac{2-2ϵ}{z(1-4z)}{I}_{14}-\frac{4-6ϵ-6(1-2ϵ)z}{z(1-4z)}{I}_{15}-\frac{6-8ϵ}{z(1-4z)}{I}_{16}+\frac{1-3ϵ}{z}{I}_{17},$$

$$I_{18}^{\prime }=-\frac{1-ϵ}{z^2(1-4z)}{I}_{14}-\frac{1-ϵ+(1-2ϵ)z}{z^2(1-4z)}{I}_{15}-\frac{3-4ϵ}{z^2(1-4z)}{I}_{16}-\frac{ϵ}{z}{I}_{18}.$$

$${\hat{G}}_{27}^{(1)2P}$$

$${\hat{G}}_{27}^{(1)3P}$$

$${\hat{G}}_{27}^{(1)m,3P}$$

$${\hat{G}}_{27}^{(1)m,2P}$$

$$\left\{f_0,f_1,g_0,g_1,r_{-1},r_0,r_1\right\}$$

$${\hat{G}}_{27}^{(1)m,3P}(z=0)=\frac{20}{27ϵ}+\frac{770}{81}+\left(\frac{18191}{243}-\frac{35}{27}{\pi }^2\right)ϵ+\mathcal{O}({ϵ}^2).$$

$$f_0(z)=-\frac{1942}{243}+2\mathrm{Re}[a(z)+b(z)].$$

$$f_1(z)$$

Full text

TTP17-007

IFT-2/2017

**NNLO QCD counterterm contributions to $\bar{B}\to X_{s}\gamma$ for the physical value of $m_{c}$

**

Mikołaj Misiak1,2, Abdur Rehman1,3 and  Matthias Steinhauser4

1 *Institute of Theoretical Physics, Faculty of Physics, University of Warsaw,

02-093 Warsaw, Poland.

2 *Theoretical Physics Department, CERN, CH-1211 Geneva 23, Switzerland.

3 *National Centre for Physics, Quaid-i-Azam University Campus,

Islamabad 45320, Pakistan.

4 *Institut für Theoretische Teilchenphysik, Karlsruhe Institute of Technology (KIT),

76128 Karlsruhe, Germany.

**Abstract

**

One of the most important ${\mathcal{O}}(\alpha_{\mathrm{s}}^{2})$ corrections to the $\bar{B}\to X_{s}\gamma$ branching ratio originates from interference of contributions from the current-current and photonic dipole operators. Its value has been estimated using an interpolation in the charm quark mass between the known results at $m_{c}=0$ and for $m_{c}\gg\frac{1}{2}m_{b}$. An explicit calculation for the physical value of $m_{c}$ is necessary to remove the associated uncertainty. In the present work, we evaluate all the ultraviolet counterterm contributions that are relevant for this purpose.

1 Introduction

The inclusive decay $\bar{B}\to X_{s}\gamma$ belongs to the most interesting flavour changing processes. Many popular extensions of the Standard Model (SM) receive important constraints from its branching ratio ${\mathcal{B}}_{s\gamma}$. For instance, the 95% C.L. bound it provides on the charged Higgs boson mass in the Two-Higgs-Doublet Model II is in the vicinity of $580\,$GeV [1, 2]. Strengthening such bounds in the future is likely to be achieved with the upcoming precise measurements of ${\mathcal{B}}_{s\gamma}$ at Belle II [3], as well as with more accurate theoretical calculations within the SM.

The SM prediction for the CP- and isospin-averaged branching ratio ${\mathcal{B}}_{s\gamma}$ amounts to [4, 5]

[TABLE]

for $E_{\gamma}>E_{0}=1.6\,$GeV. The corresponding experimental world average evaluated by HFAG [6] reads

[TABLE]

It is based on the Belle [1, 7], Babar [8, 9, 10] and CLEO [11] measurements performed at $E_{0}\in[1.7,2.0]\,$GeV, and then extrapolated down to $E_{0}=1.6\,$GeV using the method of Ref. [12]. A recent discussion of the energy extrapolation issue and its effect on the averages can be found in Ref. [2]. In that paper

[TABLE]

was obtained by naively averaging the measurements at $E_{0}=1.9\,$GeV only [1, 7, 8, 9, 10], and then performing an extrapolation to $E_{0}=1.6\,$GeV.

As far as the SM prediction is concerned, it is based on the fact that the hadronic decay rate $\Gamma(\bar{B}\to X_{s}\gamma)$ is well approximated by the corresponding partonic one $\Gamma(b\to X_{s}^{p}\gamma)$, where $X_{s}^{p}=s,sg,sgg,sq\bar{q},\ldots$. The difference between the two rates is due to non-perturbative effects. Such effects contribute at around $+3\%$ level [13] to the central value in Eq. (1.1), while the associated uncertainty is estimated [14] to be around $\pm 5\%$. This uncertainty is combined in quadrature with the remaining three ones that are due to the input parameters ($\pm 2\%$), unknown higher-order ${\mathcal{O}}(\alpha_{\mathrm{s}}^{3})$ perturbative effects ($\pm 3\%$), as well an interpolation in the charm quark mass $m_{c}$ ($\pm 3\%$). The interpolation is used to estimate some of the most important ${\mathcal{O}}(\alpha_{\mathrm{s}}^{2})$ corrections.

The purpose of our present work is providing a contribution to removing the latter uncertainty. Perturbative corrections to $\Gamma(b\to X_{s}^{p}\gamma)$ are most conveniently analyzed in the framework of an effective theory that arises after decoupling the $W$ boson and all the heavier particles. The relevant weak interaction Lagrangian is then given by a linear combination of four-quark and dipole-type operators:  ${\mathcal{L}}_{\rm weak}\sim\sum_{i}C_{i}(\mu)Q_{i}$ (see, e.g., Eq. (1.6) of Ref. [5]). For most of our discussion here, only three of them are going to be relevant, namely

[TABLE]

The $\overline{\rm MS}$-renormalized Wilson coefficients $C_{i}(\mu)$ are first evaluated at the electroweak scale $\mu_{0}\sim M_{W},m_{t}$, and then evolved down to $\mu_{b}\sim m_{b}$ using the effective theory renormalization group equations. A complete ${\mathcal{O}}(\alpha_{\mathrm{s}}^{2})$ calculation of all the relevant $C_{i}(\mu_{b})$ was finalized a decade ago, with the latest contribution amounting to a determination of the necessary four-loop anomalous dimension matrix [15].

Once the Wilson coefficients $C_{i}(\mu_{b})$ are found, one writes the perturbative rate as111 For simplicity, we ignore here the higher-order electroweak and ${\mathcal{O}}(V_{ub})$ corrections. They are included in the numerical evaluation of the SM prediction (1.1). Moreover, we make use of the fact that $C_{i}\in\mathbb{R}$ in the SM.

[TABLE]

where $G_{F}$ is the Fermi constant, and $V_{ij}$ are the Cabibbo-Kobayashi-Maskawa matrix elements. Both the $b$-quark mass $m_{b}$ and the electromagnetic coupling $\alpha_{\mathrm{e}m}$ are assumed to be on-shell renormalized.

The interference terms $\hat{G}_{ij}$ form a symmetric $8\times 8$ matrix. Their perturbative expansion in $\widetilde{\alpha}_{\mathrm{s}}\equiv\alpha_{\mathrm{s}}(\mu_{b})/(4\pi)$ reads

[TABLE]

The Next-to-Next-to-Leading Order (NNLO) corrections to be studied in the present work are given by $\hat{G}^{(2)}_{17}$ and $\hat{G}^{(2)}_{27}$. They describe interferences of decay amplitudes generated by the current-current operators $Q_{1,2}$ and the photonic dipole operator $Q_{7}$. Such interference terms are most conveniently represented via Cutkosky rules as four-loop propagator diagrams with unitarity cuts. Sample diagrams are presented in Fig. 1. Altogether, around 850 of such four-loop diagrams need to be evaluated.

Despite the fact that $\hat{G}_{(1,2)7}$ determine some of the most important NNLO corrections to ${\mathcal{B}}_{s\gamma}$ (comparable to the present experimental and theoretical uncertainties), their explicit form for the physical value of $m_{c}$ remains unknown, which is due to extreme complexity of the necessary calculation. What has been found so far for arbitrary $m_{c}$ are contributions from diagrams with massless and massive quark loops on the gluon lines [16, 17, 18]. They are the basis for evaluating $\hat{G}_{(1,2)7}$ in the Brodsky-Lepage-Mackenzie (BLM) [19] approximation. As far as the non-BLM terms are concerned, they have been calculated for $m_{c}\gg\frac{1}{2}m_{b}$ [20] and at $m_{c}=0$ [5] only. Next, their values for the physical $m_{c}$ were found using an interpolation, as described in detail in Ref. [5].

The interpolated contributions affect ${\mathcal{B}}_{s\gamma}$ by $(+5\pm 3)\%$, which implies that their effects on both the central value and the overall uncertainty are sizeable. The interpolation uncertainty could be completely removed via a calculation of $\hat{G}_{(1,2)7}$ directly at the physical value of $m_{c}$. Such a calculation involves the “bare” diagrams like the ones in Fig. 1, as well as many lower-loop diagrams with insertions of ultraviolet counterterms.

In the present paper, we provide results for all the necessary counterterm contributions in the $m_{c}\neq 0$ case. Although it is only a first step towards the complete NNLO calculation of $\hat{G}_{(1,2)7}$, it is by no means a trivial one given that two mass scales are present in the relevant three-loop propagator diagrams with unitarity cuts.222 The two mass scales are given by $m_{b}$ and $m_{c}$, while the s-quark is treated as massless throughout the paper. After generating the considered interference contributions and expressing them in terms of Master Integrals (MIs), we evaluate all the MIs using several methods. In the first approach, we derive a system of Differential Equations (DEs) for the MIs, and solve it numerically starting from boundary conditions at $m_{c}\gg\frac{1}{2}m_{b}$. The latter are found using asymptotic expansions in such a limit. The relevance of this method lies in the fact that it is basically the only feasible one for the future “bare” calculation. Thus, our current task serves as a preparation for the future enterprise. Apart from the numerical approach, we also use analytical methods for determining all the necessary MIs, either in closed forms or as expansions in $z=m_{c}^{2}/m_{b}^{2}$.

Similarly to the $m_{c}=0$ calculation of Ref. [5], we restrict ourselves to the limit $E_{0}=0$. Considering such a limit is a necessary step on the way to the $E_{0}\neq 0$ case. Results in the latter case are going to be obtained in the future by subtracting the low-$E_{\gamma}$ tail from the total rate, in the same way as it was done for $\hat{G}_{77}$ [21, 22].

Our article is organized as follows. In Section 2, we define all the relevant counterterm contributions via an explicit formula for the renormalized $\hat{G}_{(1,2)7}$. Expressions for them in terms of the MIs, as well as DEs for the MIs are given in Section 3. Our numerical results for all the considered counterterm contributions are presented in Section 4. Section 5 is devoted to presenting our final results as expansions in $z$. We conclude in Section 6.

2 The NNLO renormalization formula for $\hat{G}_{17}$ and $\hat{G}_{27}$

The renormalized $\hat{G}^{(2)}_{i7}$ for $i=1,2$ are obtained from the bare ones $\hat{G}_{i7}^{(2)\rm bare}$ and the ultraviolet counterterms. The necessary expression in the $m_{c}=0$ case together with all its ingredients was presented in Section 2.2 of Ref. [5]. Following that article, we skip the contributions that are due to charm quark loops on the gluon lines together with the corresponding counterterms. Such NNLO contributions have been already calculated for arbitrary $m_{c}$ in Ref. [18], and are included in the phenomenological analysis of Ref. [5] (Eq. (3.8) there). Under such assumptions, the renormalization formula for arbitrary $m_{c}$ takes the following form (for $i=1,2$):

[TABLE]

It turns out to differ from the $m_{c}=0$ case only in the last line where the renormalization of $m_{c}$ is taken into account. Moreover, we have simplified a few details of the notation without changing the actual content of the equation. Its l.h.s. corresponds to such a renormalization scheme where the Wilson coefficients, $\alpha_{s}$ and $m_{c}$ are $\overline{\rm MS}$-renormalized at the scale $\mu_{b}$. On the other hand, the $b$-quark mass and the external quark fields are on-shell renormalized, which is indicated by the superscript ${\rm OS}$ at the corresponding renormalization constants. In these constants (Eq. (2.8) of Ref. [5]), as well as in Eq. (2.1) here, one should substitute $s=\mu_{b}^{2}/m_{b}^{2}$. The remaining renormalization constants are the ${\rm MS}$-scheme ones (Eq. (2.9) of Ref. [5]). In particular $Z_{m}=1-4\widetilde{\alpha}_{\mathrm{s}}/\epsilon+{\mathcal{O}}(\widetilde{\alpha}_{\mathrm{s}}^{2})$  in  $D=4-2\epsilon$  spacetime dimensions.

The bare $\hat{G}_{kl}$ on the r.h.s. of Eq. (2.1) are assumed to have been evaluated for the renormalization scale $\mu^{2}=e^{\gamma}m_{b}^{2}/(4\pi)$, where $\gamma$ is the Euler-Mascheroni constant. They involve indices corresponding to the operators listed in Eqs. (1.6) and (2.5) of Ref. [5]. In some of the cases, the superscript “${\rm bare}$” has been replaced by either “$3P$” or “$m$”. In these cases, either only the three-particle cuts were included ($3P$), or one of the $b$-quark propagators was squared ($m$). The latter operation is the simplest way to obtain the counterterms that originate from renormalization of the $b$-quark mass.

Explicit expressions for all the necessary bare $\hat{G}_{kl}$ in the $m_{c}=0$ case have been given in Eqs. (2.3)-(2.7) of Ref. [5]. Only a few of them get modified in the $m_{c}\neq 0$ case, namely

[TABLE]

In the following, we shall calculate all of them except for $\hat{G}_{17}^{(2)\rm bare}$ and $\hat{G}_{27}^{(2)\rm bare}$. Some of the considered quantities are related to each other in a very simple manner. In particular

[TABLE]

It is due to the fact that in all the non-vanishing Feynman diagrams that contribute to these quantities (e.g., those in Fig. 2), the internal gluon couples only once to the charm quark loop. In consequence, the diagrams with $Q_{1}$ and $Q_{2}$ differ only by a colour factor of $-\frac{1}{6}$ that comes from the identity $T^{a}T^{b}T^{a}=-\frac{1}{6}T^{b}$ for the $SU(3)$ generators. The case is analogous for the evanescent operators

[TABLE]

Another simple identity allows us to express $\hat{G}_{7(12)}^{(1)\rm bare}$ in terms of various contributions to $\hat{G}_{27}^{(1)\rm bare}$. It reads

[TABLE]

where $\hat{G}_{27}^{(1)2P(d)}$ and $\hat{G}_{27}^{(1)2P(u)}$ stand for the two-particle cut contributions to $\hat{G}_{27}^{(1)\rm bare}$ that are proportional to the down-type and up-type quark electric charges, respectively. The identity (2.4) can be recovered after expressing all the involved quantities in terms of the MIs, but without actually calculating any of the MIs.

Thus, what remains to be calculated are the six quantities on the r.h.s. of the following two equations:

[TABLE]

Analogously to the first of the above equations, the three terms on the r.h.s. of the second one are defined according to the quark electric charges, and to the number of final-state particles.

All the considered quantities need to be found up to ${\mathcal{O}}(\epsilon)$. In the case of Eq. (2.5), it is equivalent to extending the NLO calculations of Refs. [23, 24, 25, 26, 27] to one more order in $\epsilon$. Such an extension is actually already available from Ref. [28], and we shall compare our results to that article. As far as the counterterms in Eq. (2.6) for $m_{c}\neq 0$ are concerned, our calculation is new.

3 Master integrals and differential equations

We generate the relevant Feynman diagrams using FeynArts [29]. The Dirac algebra after summing/averaging over polarizations is performed with fully anticommuting $\gamma_{5}$, which is allowed because only traces with even numbers of $\gamma_{5}$ can give non-vanishing contributions in the considered problem. Once the results are expressed in terms of scalar integrals, we reduce them to the MIs using Integration By Parts (IBP) identities, with the help of the code FIRE [30]. This way we express the six quantities on the r.h.s. of Eqs. (2.5) and (2.6) as linear combinations of eighteen master integrals $I_{k}$

[TABLE]

with coefficients $A_{k}$–$F_{k}$ that are provided in the ancillary files [31]. The MIs to be determined are displayed in Tab. 1. In our normalization conventions for them, we include an extra factor of $2^{11}\pi^{5}$ together with an appropriate integer power of $m_{b}$ that makes each $I_{k}$ dimensionless. For instance,

[TABLE]

where $p$ is the external momentum, while $k_{i}$ are the cut propagator momenta (directed from left to right in all the MIs of Tab. 1). In the above expressions, integrals over $d{\rm PS}_{n}$ stand for the $n$-body massless phase space integrals, with $n=2,3$. They can be simplified [32] to integrals over  $\hat{s}_{ij}=2(k_{i}k_{j})/m_{b}^{2}$  for arbitrary integrands $W_{2}(\hat{s}_{12})$ and $W_{3}(\hat{s}_{12},\hat{s}_{13},\hat{s}_{23})$:

[TABLE]

As already mentioned, we set $\mu^{2}=e^{\gamma}m_{b}^{2}/(4\pi)$ when evaluating the MIs. With the normalization as defined above, all the master integrals $I_{k}$ and the coefficients $A_{k},\ldots,F_{k}$ in Eq. (3.1) are functions of two variables only: $\epsilon$ and $z=m_{c}^{2}/m_{b}^{2}$.

In four spacetime dimensions, the global normalization of all the considered contributions to the decay rate is fixed by Eq. (1.5). In $D=4-2\epsilon$ dimensions, we follow Ref. [5] where the global normalization is fixed by

[TABLE]

The factor $(1-\epsilon)$ above comes from the Dirac algebra in $\left|\langle s\gamma|Q_{7}|b\rangle\right|^{2}$ summed over polarizations.

To evaluate the MIs, we differentiate them with respect to $z$, and reduce the resulting integrals to the MIs again, using the IBP identities. This gives us a closed set of DEs which takes the following form ($I^{\prime}_{k}=dI_{k}/dz$):

[TABLE]

The necessary boundary conditions for the above DEs are evaluated in the $z\gg 1$ limit, in which the integrals become considerably simpler. Large-$z$ expansions of the three-body MIs ($I_{14}$–$I_{18}$) are obtained in a straightforward manner using the Feynman parameterization. In the two-body case ($I_{1}$–$I_{13}$), we have used the code exp [33, 34] that applies the asymptotic expansion method in a fully automatic manner.

4 Numerical results for the counterterm contributions

In our first approach, we expanded all the MIs in $\epsilon$ up to the relevant orders, which gave us a system of coupled ordinary DEs for 76 functions of $z$ being a single variable. The boundary condition was chosen at $z=20$, and evaluated using the expansions up to ${\mathcal{O}}(1/z^{8})$. Since the coefficients of the DEs (3.6) have singularities at $z\in\{0,\frac{1}{4},1\}$, the variable $z$ must be treated as complex, and a contour in the complex $z$-plane must be chosen for the numerical solution. In practice, we solved the DEs along ellipses in the lower half-plane of $z$, following the Feynman prescription $m_{c}^{2}\to m_{c}^{2}-i\varepsilon$ in the propagators, and treating $m_{b}^{2}$ in $z=m_{c}^{2}/m_{b}^{2}$ as a real rescaling factor. In this way, non-vanishing imaginary parts in some of the MIs were obtained on the real axis for $z<\frac{1}{4}$ (below the $c\bar{c}$ production threshold), despite the fact that both the boundary conditions at $z=20$ and the DE coefficients were purely real. In particular, the proper imaginary parts of the functions $a(z)$ and $b(z)$ from Eqs. (3.3)–(3.4) of Ref. [27] were recovered. Such imaginary parts drop out from our final results in Eqs. (2.5)–(2.6) but, nevertheless, they provide convenient cross-checks of the calculation.

The numerical solution was obtained using the FORTRAN code ZVODE [35], upgraded to quadrupole-double precision with the help of the QD computation package [36]. The physical value of $z$ is around $0.06$. However, apart from the vicinity of this point, we have considered a wide range of the final values, in order to test the limits at $z\to 0$, as well as the behaviour around the threshold at $z=\frac{1}{4}$.

To present our numerical results for the quantities listed in Eq. (3.1), we parameterize them in terms of nine functions of $z$ as follows

[TABLE]

Their dependence on $z$ is shown in Fig. 3. Each of the (blue) dots corresponds to a particular final value of the numerical solution of the DEs along an ellipse in the complex plane. The physical point around $z=0.06$ is marked by a slightly bigger (red) dot. Similar (red) dots on the vertical axes are used to indicate the limits at $z\to 0$ in all the cases when they are finite, i.e. for all the functions except those in Eq. (4.3). They can be taken over from the calculation performed at $z=0$ from the outset [5]:

[TABLE]

where $\zeta_{k}\equiv\zeta(k)$. In the case of $\hat{G}_{27}^{(1)m,3P}$, one finds [5]

[TABLE]

The above expression contains a $1/\epsilon$ divergence, contrary to Eq. (4.3). It arises due to the logarithmic divergences of the functions $j_{i}(z)$ at $z\to 0$. Such non-commutation of limits is frequently encountered in the framework of dimensional regularization.

Our expansions above $z=20$ that have served as the boundary conditions for the DEs are displayed as solid (blue) lines in the first six plots of Fig. 3. The remaining (green) solid lines describe expansions either in $z$ (for $z<\frac{1}{4}$) or in $\frac{1}{z}$ (for $z>\frac{1}{4}$) that have been obtained in our second approach – see Section 5. For $f_{i}(z)$, the expansion plots are terminated away from the threshold at $z=\frac{1}{4}$, as they become inaccurate in its vicinity. In the three-body cases ($g_{i}(z)$, $j_{i}(z)$) the applied expansions are so accurate that their mismatch at $z=\frac{1}{4}$ is invisible within the plot resolution.333 The solid lines for $g_{0}(z)$ actually correspond to the fully analytical result of Eq. (5.3).

We refrain from presenting here any numerical fits for the functions plotted in Fig. 3. Instead, either expansions or exact expressions for them will be given in Section 5. It is evident from the plots that our expansions around $z=0$ are sufficiently precise for phenomenological purposes in the vicinity of the physical point.

5 Expansions in $z$

In this section, we collect our final results for the small-$z$ expansions of the functions defined in Eqs. (4.1)–(4.4). They have been obtained in Ref. [37] with the help of Feynman parameterizations, Mellin-Barnes techniques [38], as well as our differential equations in Eq. (3.6) that allowed to extend the expansions to higher orders.

For $f_{0}(z)$, we properly recover the NLO functions $a(z)$ and $b(z)$ known from the previous calculations [24, 25, 26, 27], namely

[TABLE]

The expansions of $a(z)$ and $b(z)$ up to ${\mathcal{O}}(z^{6})$ are given in Eqs. (3.8)–(3.9) of Ref. [27], and we shall not repeat them here. For $f_{1}(z)$, we obtain

[TABLE]

where $L=\ln z$. The above result is in perfect agreement with Eqs. (15)–(17) of Ref. [28]. As far as $g_{0}(z)$ is concerned, a fully analytical expression is available

[TABLE]

where  $s=\sqrt{1-4z}$,  $X=\ln\frac{1+s}{2}-\frac{1}{2}L$,  $t=\sqrt{4z-1}$,  and  $Y=\arctan(1/t)$. For $g_{1}(z)$, we find

[TABLE]

It cannot be derived from Eq. (23) in Ref. [28] that has been obtained with the help of phase-space integrations in four dimensions (as stated above Eq. (19) there). Our result corresponds to $D=4-2\epsilon$ throughout the calculation, as required by the so-called reverse unitarity method [39] that we have applied for the IBP in the three-body case.

The remaining five functions correspond to quantities that have not been considered in Ref. [28] because they are unrelated to the charm quark mass renormalization. Our results for their expansions in $z$ read

[TABLE]

In the latter equation, the constant $C\simeq-292.228$ has been determined only numerically.

Eqs. (5.4)–(5.5) do not explicitly show all the terms that have been taken into account for the plots in Fig. 3. In the latter case, we used the series terminating at $z^{11}$, $z^{10}$, $z^{10}$, $z^{5}$ and $z^{5}$ for $g_{1}$, $j_{0}$, $j_{1}$, $r_{0}$ and $r_{1}$, respectively. Extended small-$z$ expansions, large-$z$ expansions, as well as fully analytical results for many of our master integrals can be found in Ref. [37].

It is likely that fully analytical results for all the considered counterterm contributions could be found using the so-called canonical basis [40, 41] of the MIs, and solving the DEs order-by-order in $\epsilon$, in terms of iterated integrals. However, at the time the main part of our calculation was performed, no automatic tools for constructing canonical bases at three loops (see, e.g., [42, 43]) were available. On the other hand, the expansions in $z$ are certainly sufficient for all the phenomenological purposes in our case, while the numerical method is most suitable for preparing the bare calculation.

6 Summary

We calculated all the relevant counterterm contributions to $\hat{G}^{(2)}_{17}$ and $\hat{G}^{(2)}_{27}$ at the NNLO for the physical value of the charm quark mass $m_{c}$. In combination with the future bare calculation of these quantities, the current results are necessary for removing one of the most important perturbative uncertainties in the SM prediction ${\mathcal{B}}_{s\gamma}$. In our first approach, we derived differential equations for the master integrals, and solved them numerically using the boundary conditions at large $m_{c}$. In the second approach, a combination of various techniques was used to determine the master integrals in terms of power-logarithmic expansions around $m_{c}=0$. In most of the cases, logarithmic divergences at $m_{c}\to 0$ are absent, and our results properly converge to the $m_{c}=0$ ones previously found in Ref. [5].

Acknowledgments

We are grateful to Michał Czakon and Paul Fiedler for their advice concerning numerical solutions of the differential equations, and for sharing their codes with us. M.M. and A.R. acknowledge support from the National Science Centre (Poland) research project, decision no DEC-2014/13/B/ST2/03969. The research of M.S. has been supported by the BMBF grant 05H15VKCCA.