Two loop virtual corrections to $b\to(d,s)\ell^+\ell^-$ and $c\to u\ell^+\ell^-$ for arbitrary momentum transfer

Non-factorizable two loop corrections to heavy to light flavor changing neutral current transitions due to matrix elements of current-current operators are calculated analytically for arbitrary momentum transfer. This extends previous works on $b\to(d,s)\ell^+\ell^-$ transitions. New results for $c\to u\ell^+\ell^-$ transitions are presented. Recent works on polylogarithms are used for the master integrals. For $b\to s\ell^+\ell^-$ transitions, the corrections are most significant for the imaginary parts of the effective Wilson coefficients in the large hadronic recoil range. Analytical results and ready-to-use fitted results for a specific set of parameters are provided.


Introduction
Recently, several discrepancies with the standard model (SM) have been revealed in b → s + − induced decays, e.g., [1]. However, those that do not involve lepton flavor universality violation are not definitely settled due to poorly controlled non-perturbative effects of quantum chromodynamics (QCD). Exemplary, non-factorizable corrections to form factors are commonly considered as one unreducible uncertainty while interpreting beyond the standard model (BSM) physics. Furthermore, in c → u + − transitions, nonfactorizable corrections yield the largest perturbative contribution due to the Glashow-Iliopolus-Maiani (GIM) cancellation of factorizable contributions [2][3][4]. In this paper, we improve on the current state by calculating the two loop virtual corrections to any heavy to light transition, including c and b decays, induced by current-current operators for arbitrary invariant dilepton mass q 2 .
Partonic transitions are a first approximation to the corresponding inclusive hadronic decays in the framework of an operator product expansion (OPE). This approximation is applicable away from resonance regions and up to power corrections. The resonance regions can be handled by appropriate kinematical cuts and power corrections can be treated within a heavy quark expansion (an expansion in inverse powers of the heavy quark mass, see [5] for b → s transitions). Furthermore, perturbative results are the basis for estimates of non-perturbative effects in inclusive and exclusive hadronic decays, e.g., [6,7].
Several calculations were performed for b → (d, s) transitions: In [8] and [9] b → s and b → d transitions, respectively, were computed for small q 2 . The calculation of b → s transitions for large q 2 was accomplished in [10]. A seminumerical approach was employed in [5] to present results based on b → s transitions for the full q 2 range. In [11] b → d transitions for any q 2 range were computed, extending [9]. For c → u transitions, see [12]. We emphasize that available results for b → (d, s) transitions only cover different limits, i.e. zero masses, small and large q 2 ranges, and that results for c → u transitions have become available only recently [4,12]. The analytic calculation presented in this paper covers the full q 2 range, arbitrary masses and electric charges.
Generally, the effective weak Lagrangian for heavy to light quark flavor changing neutral current (FCNC) transitions h → l at the low-energy scale µ is written as where the sum is over light quark fields q with masses below µ, and products of Cabibbo- and c → u transitions, respectively. Here C i are the Wilson coefficients and the physical operators P i , which are relevant for this paper, read where , and T a are the SU (3) C generators normalized to Tr[T a T b ] = δ ab 2 . Furthermore, G F is the Fermi constant, F µ 1 µ 2 denotes the electromagnetic field strength tensor, and g s and e are the strong and electromagnetic coupling constants, respectively.
We calculate the two loop QCD matrix elements of P 1/2 , in terms of form factors (i.e., for inclusive decays, effective Wilson coefficients) multiplying the matrix elements of P 7,9 , for h → l + − transitions. The result is valid for a general class of heavy to light transitions with arbitrary invariant momentum transfer and masses, when the mass of the light quark is neglected. This includes the transitions b → (d, s) via (uū, cc) loops, and c → u via (dd, ss) loops, where the loop quark-antiquark pair (qq) is annihilated and a photon is emitted, which may then couple to a lepton pair. Note that the computation of two loop matrix elements presented in this paper is not restricted to SM applications, see [12] for an example in leptoquark models. We use the recent works [13] and [14] for the master integrals (MIs) and their numerical evaluation, respectively.
We outline our calculation in the section 2, see also [12]. The numerical evaluation is detailed in section 2.1. Results are given in section 3, where we also comment on the phenomenological impact for b → (d, s) transitions. For the phenomenology of c → u transitions, see [12]. Appendix A contains a description of supplemented files, which encode the results of this paper.

Outline of the calculation
In this section, we outline the calculation of the diagrams shown in figure 1. Each of the five subsets represents a gauge invariant class. A sixth class, shown in figure 2, preserves the operator structure of P 9 , hence it is commonly considered as a correction to the matrix element of P 9 [8] and not included in the calculation presented in this paper. Nevertheless, we have checked that in this class only the diagram with a photon coupling to the loop of the quark-antiquark pair is non-zero and factorizes into two one loop integrals. It is the only diagram with infrared and collinear singularities.
We calculate the diagrams in figure 1 with insertions of P 2 . Insertions of P 1 are then given by color factors due to additional generators in the definition of the operator: The expressions for the first four and last (two) classes are multiplied with −1 6 and 4 3 , respectively. The matrix element for an off-shell photon γ * can be decomposed as [11] where q µ = (p h − p l ) µ is the transferred momentum, and X µ is the sum of the amputated diagrams in figure 1.
The scalar form factors F (i) are given as [11] with the projectors and coefficients C ij . The factors ( / p h + m h ) and / p l reflect the on-shell conditions for the external quarks. The form factors for l |P 2 |h , thus the effective Wilson coefficients, are proportional to F (7,9) . By gauge invariance, the form factor F (q) vanishes for each class, which we have checked.
For the calculation we utilize the computer programs FORM 4.0 [15] and REDUZE 2 [16]. We use the former for algebraic manipulations, e.g. the tensor algebra, and the latter is used for the reduction to MIs. The program REDUZE 2 is based on the Laporta algorithm [17] which employs integration by parts (IBP) identities [18] and Lorentz invariance (LI) identities [19]. Indeed, LI identities can be written as linear combinations of IBP identities [20]. We have calculated relations for each diagram based on LI identities by hand and checked them against the reduction tables built with REDUZE 2.
We find that the following diagrams in figure 1 vanish: In the fifth class, both diagrams with photons attached to the external lines vanish. The diagram with a photon emitted from the heavy quark line in the first class is zero. For F (7) all diagrams with photons attached to the external lines vanish. Furthermore, all diagrams of the fifth class vanish for F (7) . This implies that the P 1 induced dipole form factor is given by As for the set of the MIs, we match a subset of the integrals calculated in [13]. A canonical set is given in [10], where the MIs are calculated for large q 2 , yet the set of integrals is not minimal. While matching the set of [10] onto the MIs taken from [13], we find additional relations among the former, e.g. for the integral I d23 and the one of equation (A12) in [10]. Furthermore, we do not encounter the integrals of equations (A5), (A8), and (A14) in [10].
With the MIs from [13] the unrenormalized form factors are expressed in terms of (generalized) harmonic polylogarithms (HPLs), see the next section for the numerical evaluation. The prescription for the renormalization is described in, e.g., [8]. Accordingly, the operator renormalization constants are written as where the dimension is 4 − 2 . Extending the set of the physical operators P i by the evanescent operators E 11,12 defined as the coefficients a ij are compactly written as and Here P 4 = (l L γ µ 1 T a h L ) {q:mq<µ} (qγ µ 1 T a q), n f is the number of flavors and q e,i are the charges of the external and internal quarks, respectively, i.e. q e = −1 3 , q i = 2 3 for b → (d, s) and q e = 2 3 , q i = −1 3 for c → u transitions. The coefficients a 11,22 ij can be obtained from the leading order (LO) anomalous dimension matrix (ADM) γ (0) and the coefficients a 12 ij from [21], where γ (1) is the next-to leading order (NLO) ADM, and the mixing via the evanescent operators is found aŝ The counterterm form factors F ct (7,9) i→ (7,9) , due to the mixing of P 1/2 into P 7,9 , and the one loop renormalization of g s in the definition of P 9 , are [8,9] F ct(7) where β 0 = 11 − 2 3 n f . The counterterms F ct (7,9) i→4quark , due to the mixing of P 1/2 into four-quark operators P j , are defined by [8] We calculate them to O( ) and for insertions of P 1,2,4 , E 11,12 according to eqs. (18) and (12), respectively. We write where in the massless limit, m 2 q → 0, and the residual m 2 q dependence cancels to O( ) in eq. (19). With this, the matrix elements are evaluated to We have checked that expanding the matrix elements in small q 2 and in the limit m 2 q → 0 yields the results in [8] and [11], respectively.
Following [8], the renormalization of the mass m q is given by the replacement m q → Z mq m q in l |P 1,2 |h 1loop of eq. (21). The mass renormalization constants Z m in the modified minimal subtraction (MS) and the pole mass scheme are [22] Expanding the matrix elements at O( αs 4π ) and O( 0 ) gives the counterterms F ct (9) i,mq and F ct (7) i,mq = 0. We have checked for both schemes that expanding the counterterm F (9) i,mq in small q 2 yields the results in [8].
Finally, the renormalized form factors are given by subtracting F ct(7,9) i = (F ct (7,9) i→7,9 + F ct (7,9) i→4quark + F ct (7,9) i,mq ) from the unrenormalized form factors. Note that wave function renormalization of the external quark fields would need to be taken into account only if we wanted to include the diagrams of figure 2 [8]. We have checked that F ct(7) i | 0 agrees with the results in [5].

Numerical evaluation
The analytical results for the renormalized form factors as found in the previous section are provided as supplemented files, see appendix A. In this section, we describe the numerical evaluation of the MIs [13], which are expressed in terms of HPLs. Note that the undetermined M19tildeep4 [13] in the MIs cancels in the form factors. In accordance with [13], we choose the analytical continuation by subtracting from the internal mass an infinitesimal imaginary part η > 0, i.e. m 2 i → m 2 i − iη. We write the HPLs as generalized/Goncharov polylogarithms (GPLs), see, e.g., [13] (and references therein) and feed them into the computer package lieevaluate [14]. Other packages for the numerical evaluation of GPLs can be found in [23,24].
While evaluating the expressions with lieevaluate, we encounter the undefined function θ(0) and the function MyP [14] that is divergent for some arguments of GPLs [14]. We regulate the former by a perturbation of these arguments, see [14]. We have numerically checked that the form factors are insensitive to the actual choice of such perturbations. The divergences due to the function MyP cancel in the form factors within the numerical precision.
We have checked that the numerical evaluation of the MIs with lieevaluate yields a precision better than 10 −6 with respect to [13]. 1 Moreover, numerical agreement of the MIs in terms of HPLs is found by use of the package HypExp 2 [25,26] and also via numerical integration within Mathematica, yet the convergence is partially slow.
The analytical expressions are lengthy and their numerical evaluation is involved. Hence, we evaluate the form factors for fixed mass parameters and at different q 2 points.
Note that the numerical evaluation close to the kinematical endpoints q 2 = 0, m 2 h is sensitive to the ratio of mass parameters and q 2 . Finally, we fit the points.
Subtracting the counterterm form factors from the unrenormalized form factors, the 1 2 and 1 divergences cancel numerically. We have checked our calculation against the ones of [8,10,11] for b → (d, s) transitions, finding numerical agreement for different q 2 , scales, mass schemes, and parameters, see the next section. For c → uγ transitions, the effective dipole Wilson coefficient at q 2 = 0 induced by P 2 was calculated in [2]. Adding the constant terms given in [3] we have checked the calculations, finding numerical agreement.

Results
Our fitted results of the renormalized form factors F m u,d ≈ 0 and set µ = m h if not stated otherwise. We note the following: • The results of [8,10,11] agree well with our results, hence the former are only partially visible in the plots.
• The form factors are divergent at the internal quark pair mass squared q 2 = (2m i ) 2 .
• The lower left plot of figure 3 indicates a numerical precision of our results better than 10 −3 with respect to the analytical result of [11].
• The lower right plot of figure 3 shows agreement with the expansion of [10] at high q 2 . Note that we do not plot the ratio close to the kinematical endpoint q 2 = m 2 b , where the evaluation of the result of [10] yields oscillations. At low q 2 , a breakdown of the convergence of the expansion of [8] is visible around q 2 ∼ 5 GeV 2 , where the sensitivity on µ increases. Similar conclusions are drawn from figure 5.
Furthermore, note that our definition of the form factors and the one in [8,10] differs by a minus sign. Recall that F 2 . For F 1,2 , the numerical precision of our results is better than 10 −2 with respect to the results of [8,10,11] for b → (d, s) transitions.
We deduce that the overall numerical precision of our result is better than a percent. To have ready-to-use result at hand, we provide our results for the masses of eq. (23) as supplemented files, see appendix A.
Next, we comment on the impact of the new results for b → (d, s) transitions. For massless internal quarks an independent analytical result is provided in [11]. The high q 2 range for massive internal quarks is well approximated by the results of [10]. On the other side, compared to the low q 2 range [8] for massive internal quarks, our results indicate significant corrections around q 2 ∼ 5 GeV 2 . Following [5] and [10] for the matrix element     [8,10], whereas, in the left plots, the results of [11] agree with the solid lines.  Re[C eff ] [8] (blue) and Im[C eff ] this Im[C eff ] [8] (purple) of the effective Wilson coefficients C eff 7 (left) andC eff 9 (right) for b → (d, s) transitions induced by a massive internal charm quark. We follow [5,10], adding the results of this paper and the expansion of [8] for µ = 5 GeV. A ratio of one is indicated by the black dotted line. of the chromomagnetic operator, the effective Wilson coefficientsC eff 7,9 at next-to-next-to leading logarithmic (NNLL) order are shown in figures 7, 8. We add our results and, for comparison, the results of [8] for 1 GeV 2 ≤ q 2 ≤ 7 GeV 2 . Note that for phenomenological purposesP 7,9 = αs 4π P 7,9 , thusC 7,9 = 4π αs C 7,9 . One observes that the real parts are stable, whereas corrections to imaginary parts increase to several percent towards higher q 2 . The effective Wilson coefficients obey the hierarchies Re[C eff 7,9 ] Im[C eff 7,9 ]. Thus, observables which only depend on the real parts or magnitudes of the effective Wilson coefficients are marginally affected by our results as long as the results of [8,10,11] are taken into account. On the other hand, effects on observables which depend on the imaginary parts of the effective Wilson coefficients are significant in the low q 2 range.

Summary
In this paper, we calculated the effective Wilson coefficients for heavy to light FCNC transitions induced by the matrix elements of current-current operators at two loop. The results are valid for arbitrary momentum transfer and masses, when the light mass is neglected. For the MIs we used the works of [13,14]. Our computation extends previous works on b → (d, s) + − transitions [8][9][10][11] and is new for c → u + − transitions. We found significant corrections to the imaginary parts of the effective Wilson coefficients for b → (d, s) transitions in the large hadronic recoil range, see figure 8, which should be included in future analyses, including BSM studies. Other corrections for b → (d, s) transitions were found to be marginal. For c → u transitions and an application to leptoquark models, we refer to [12]. Along with this paper we provide supplemented files, containing our analytical results and fitted results for a specific set of parameters. Finally, our calculation is an independent check of the results of [13,14] and [8][9][10][11], with which we agree in the corresponding limits.